LCOV - code coverage report
Current view: top level - src/aobasis - basis_set_types.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:e7e05ae) Lines: 1360 1441 94.4 %
Date: 2024-04-18 06:59:28 Functions: 27 32 84.4 %

          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             : !> \par History
      10             : !>      - 02.2004 flexible normalization of basis sets [jgh]
      11             : !>      - 07.2014 Add a set of contraction coefficient that only work on active
      12             : !>                functions
      13             : !> \author Matthias Krack (04.07.2000)
      14             : ! **************************************************************************************************
      15             : MODULE basis_set_types
      16             : 
      17             :    USE ai_coulomb,                      ONLY: coulomb2
      18             :    USE bibliography,                    ONLY: VandeVondele2007,&
      19             :                                               cite_reference
      20             :    USE cp_linked_list_input,            ONLY: cp_sll_val_next,&
      21             :                                               cp_sll_val_type
      22             :    USE cp_parser_methods,               ONLY: parser_get_object,&
      23             :                                               parser_search_string
      24             :    USE cp_parser_types,                 ONLY: cp_parser_type,&
      25             :                                               parser_create,&
      26             :                                               parser_release
      27             :    USE input_section_types,             ONLY: section_vals_list_get,&
      28             :                                               section_vals_type,&
      29             :                                               section_vals_val_get
      30             :    USE input_val_types,                 ONLY: val_get,&
      31             :                                               val_type
      32             :    USE kinds,                           ONLY: default_path_length,&
      33             :                                               default_string_length,&
      34             :                                               dp
      35             :    USE mathconstants,                   ONLY: dfac,&
      36             :                                               pi
      37             :    USE memory_utilities,                ONLY: reallocate
      38             :    USE message_passing,                 ONLY: mp_para_env_type
      39             :    USE orbital_pointers,                ONLY: coset,&
      40             :                                               indco,&
      41             :                                               init_orbital_pointers,&
      42             :                                               nco,&
      43             :                                               ncoset,&
      44             :                                               nso,&
      45             :                                               nsoset
      46             :    USE orbital_symbols,                 ONLY: cgf_symbol,&
      47             :                                               sgf_symbol
      48             :    USE orbital_transformation_matrices, ONLY: orbtramat
      49             :    USE sto_ng,                          ONLY: get_sto_ng
      50             :    USE string_utilities,                ONLY: integer_to_string,&
      51             :                                               remove_word,&
      52             :                                               uppercase
      53             :    USE util,                            ONLY: sort
      54             : #include "../base/base_uses.f90"
      55             : 
      56             :    IMPLICIT NONE
      57             : 
      58             :    PRIVATE
      59             : 
      60             :    ! Global parameters (only in this module)
      61             : 
      62             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'basis_set_types'
      63             : 
      64             :    ! basis set sort criteria
      65             :    INTEGER, PARAMETER, PUBLIC               :: basis_sort_default = 0, &
      66             :                                                basis_sort_zet = 1
      67             : 
      68             : ! **************************************************************************************************
      69             :    ! Define the Gaussian-type orbital basis set type
      70             : 
      71             :    TYPE gto_basis_set_type
      72             :       !MK PRIVATE
      73             :       CHARACTER(LEN=default_string_length)       :: name = ""
      74             :       CHARACTER(LEN=default_string_length)       :: aliases = ""
      75             :       REAL(KIND=dp)                              :: kind_radius = 0.0_dp
      76             :       REAL(KIND=dp)                              :: short_kind_radius = 0.0_dp
      77             :       INTEGER                                    :: norm_type = -1
      78             :       INTEGER                                    :: ncgf = -1, nset = -1, nsgf = -1
      79             :       CHARACTER(LEN=12), DIMENSION(:), POINTER   :: cgf_symbol => NULL()
      80             :       CHARACTER(LEN=6), DIMENSION(:), POINTER    :: sgf_symbol => NULL()
      81             :       REAL(KIND=dp), DIMENSION(:), POINTER       :: norm_cgf => NULL(), set_radius => NULL()
      82             :       INTEGER, DIMENSION(:), POINTER             :: lmax => NULL(), lmin => NULL(), &
      83             :                                                     lx => NULL(), ly => NULL(), lz => NULL(), &
      84             :                                                     m => NULL(), ncgf_set => NULL(), &
      85             :                                                     npgf => NULL(), nsgf_set => NULL(), nshell => NULL()
      86             :       REAL(KIND=dp), DIMENSION(:, :), POINTER    :: cphi => NULL(), pgf_radius => NULL(), sphi => NULL(), &
      87             :                                                     scon => NULL(), zet => NULL()
      88             :       INTEGER, DIMENSION(:, :), POINTER          :: first_cgf => NULL(), first_sgf => NULL(), l => NULL(), &
      89             :                                                     last_cgf => NULL(), last_sgf => NULL(), n => NULL()
      90             :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc => NULL()
      91             :    END TYPE gto_basis_set_type
      92             : 
      93             :    TYPE gto_basis_set_p_type
      94             :       TYPE(gto_basis_set_type), POINTER :: gto_basis_set => NULL()
      95             :    END TYPE gto_basis_set_p_type
      96             : 
      97             : ! **************************************************************************************************
      98             :    ! Define the Slater-type orbital basis set type
      99             : 
     100             :    TYPE sto_basis_set_type
     101             :       PRIVATE
     102             :       CHARACTER(LEN=default_string_length)       :: name = ""
     103             :       INTEGER                                    :: nshell = -1
     104             :       CHARACTER(LEN=6), DIMENSION(:), POINTER    :: symbol => NULL()
     105             :       INTEGER, DIMENSION(:), POINTER             :: nq => NULL(), lq => NULL()
     106             :       REAL(KIND=dp), DIMENSION(:), POINTER       :: zet => NULL()
     107             :    END TYPE sto_basis_set_type
     108             : 
     109             : ! **************************************************************************************************
     110             :    INTERFACE read_gto_basis_set
     111             :       MODULE PROCEDURE read_gto_basis_set1, read_gto_basis_set2
     112             :    END INTERFACE
     113             : ! **************************************************************************************************
     114             : 
     115             :    ! Public subroutines
     116             :    PUBLIC :: allocate_gto_basis_set, &
     117             :              deallocate_gto_basis_set, &
     118             :              get_gto_basis_set, &
     119             :              init_aux_basis_set, &
     120             :              init_cphi_and_sphi, &
     121             :              init_orb_basis_set, &
     122             :              read_gto_basis_set, &
     123             :              copy_gto_basis_set, &
     124             :              create_primitive_basis_set, &
     125             :              combine_basis_sets, &
     126             :              set_gto_basis_set, &
     127             :              sort_gto_basis_set, &
     128             :              write_gto_basis_set, &
     129             :              write_orb_basis_set
     130             : 
     131             :    PUBLIC :: allocate_sto_basis_set, &
     132             :              read_sto_basis_set, &
     133             :              create_gto_from_sto_basis, &
     134             :              deallocate_sto_basis_set, &
     135             :              set_sto_basis_set, &
     136             :              srules
     137             : 
     138             :    ! Public data types
     139             :    PUBLIC :: gto_basis_set_p_type, &
     140             :              gto_basis_set_type, &
     141             :              sto_basis_set_type
     142             : 
     143             : CONTAINS
     144             : 
     145             : ! **************************************************************************************************
     146             : !> \brief ...
     147             : !> \param gto_basis_set ...
     148             : ! **************************************************************************************************
     149       19164 :    SUBROUTINE allocate_gto_basis_set(gto_basis_set)
     150             : 
     151             :       ! Allocate a Gaussian-type orbital (GTO) basis set data set.
     152             : 
     153             :       ! - Creation (26.10.2000,MK)
     154             : 
     155             :       TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set
     156             : 
     157       19164 :       CALL deallocate_gto_basis_set(gto_basis_set)
     158             : 
     159       19164 :       ALLOCATE (gto_basis_set)
     160             : 
     161       19164 :    END SUBROUTINE allocate_gto_basis_set
     162             : 
     163             : ! **************************************************************************************************
     164             : !> \brief ...
     165             : !> \param gto_basis_set ...
     166             : ! **************************************************************************************************
     167       38620 :    SUBROUTINE deallocate_gto_basis_set(gto_basis_set)
     168             : 
     169             :       ! Deallocate a Gaussian-type orbital (GTO) basis set data set.
     170             : 
     171             :       ! - Creation (03.11.2000,MK)
     172             : 
     173             :       TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set
     174             : 
     175       38620 :       IF (ASSOCIATED(gto_basis_set)) THEN
     176       19456 :          IF (ASSOCIATED(gto_basis_set%cgf_symbol)) DEALLOCATE (gto_basis_set%cgf_symbol)
     177       19456 :          IF (ASSOCIATED(gto_basis_set%sgf_symbol)) DEALLOCATE (gto_basis_set%sgf_symbol)
     178       19456 :          IF (ASSOCIATED(gto_basis_set%norm_cgf)) DEALLOCATE (gto_basis_set%norm_cgf)
     179       19456 :          IF (ASSOCIATED(gto_basis_set%set_radius)) DEALLOCATE (gto_basis_set%set_radius)
     180       19456 :          IF (ASSOCIATED(gto_basis_set%lmax)) DEALLOCATE (gto_basis_set%lmax)
     181       19456 :          IF (ASSOCIATED(gto_basis_set%lmin)) DEALLOCATE (gto_basis_set%lmin)
     182       19456 :          IF (ASSOCIATED(gto_basis_set%lx)) DEALLOCATE (gto_basis_set%lx)
     183       19456 :          IF (ASSOCIATED(gto_basis_set%ly)) DEALLOCATE (gto_basis_set%ly)
     184       19456 :          IF (ASSOCIATED(gto_basis_set%lz)) DEALLOCATE (gto_basis_set%lz)
     185       19456 :          IF (ASSOCIATED(gto_basis_set%m)) DEALLOCATE (gto_basis_set%m)
     186       19456 :          IF (ASSOCIATED(gto_basis_set%ncgf_set)) DEALLOCATE (gto_basis_set%ncgf_set)
     187       19456 :          IF (ASSOCIATED(gto_basis_set%npgf)) DEALLOCATE (gto_basis_set%npgf)
     188       19456 :          IF (ASSOCIATED(gto_basis_set%nsgf_set)) DEALLOCATE (gto_basis_set%nsgf_set)
     189       19456 :          IF (ASSOCIATED(gto_basis_set%nshell)) DEALLOCATE (gto_basis_set%nshell)
     190       19456 :          IF (ASSOCIATED(gto_basis_set%cphi)) DEALLOCATE (gto_basis_set%cphi)
     191       19456 :          IF (ASSOCIATED(gto_basis_set%pgf_radius)) DEALLOCATE (gto_basis_set%pgf_radius)
     192       19456 :          IF (ASSOCIATED(gto_basis_set%sphi)) DEALLOCATE (gto_basis_set%sphi)
     193       19456 :          IF (ASSOCIATED(gto_basis_set%scon)) DEALLOCATE (gto_basis_set%scon)
     194       19456 :          IF (ASSOCIATED(gto_basis_set%zet)) DEALLOCATE (gto_basis_set%zet)
     195       19456 :          IF (ASSOCIATED(gto_basis_set%first_cgf)) DEALLOCATE (gto_basis_set%first_cgf)
     196       19456 :          IF (ASSOCIATED(gto_basis_set%first_sgf)) DEALLOCATE (gto_basis_set%first_sgf)
     197       19456 :          IF (ASSOCIATED(gto_basis_set%l)) DEALLOCATE (gto_basis_set%l)
     198       19456 :          IF (ASSOCIATED(gto_basis_set%last_cgf)) DEALLOCATE (gto_basis_set%last_cgf)
     199       19456 :          IF (ASSOCIATED(gto_basis_set%last_sgf)) DEALLOCATE (gto_basis_set%last_sgf)
     200       19456 :          IF (ASSOCIATED(gto_basis_set%n)) DEALLOCATE (gto_basis_set%n)
     201       19456 :          IF (ASSOCIATED(gto_basis_set%gcc)) DEALLOCATE (gto_basis_set%gcc)
     202       19456 :          DEALLOCATE (gto_basis_set)
     203             :       END IF
     204       38620 :    END SUBROUTINE deallocate_gto_basis_set
     205             : 
     206             : ! **************************************************************************************************
     207             : !> \brief ...
     208             : !> \param basis_set_in ...
     209             : !> \param basis_set_out ...
     210             : ! **************************************************************************************************
     211        2940 :    SUBROUTINE copy_gto_basis_set(basis_set_in, basis_set_out)
     212             : 
     213             :       ! Copy a Gaussian-type orbital (GTO) basis set data set.
     214             : 
     215             :       TYPE(gto_basis_set_type), INTENT(IN)               :: basis_set_in
     216             :       TYPE(gto_basis_set_type), POINTER                  :: basis_set_out
     217             : 
     218             :       INTEGER                                            :: maxco, maxpgf, maxshell, ncgf, nset, nsgf
     219             : 
     220        2940 :       CALL allocate_gto_basis_set(basis_set_out)
     221             : 
     222        2940 :       basis_set_out%name = basis_set_in%name
     223        2940 :       basis_set_out%aliases = basis_set_in%aliases
     224        2940 :       basis_set_out%kind_radius = basis_set_in%kind_radius
     225        2940 :       basis_set_out%norm_type = basis_set_in%norm_type
     226        2940 :       basis_set_out%nset = basis_set_in%nset
     227        2940 :       basis_set_out%ncgf = basis_set_in%ncgf
     228        2940 :       basis_set_out%nsgf = basis_set_in%nsgf
     229        2940 :       nset = basis_set_in%nset
     230        2940 :       ncgf = basis_set_in%ncgf
     231        2940 :       nsgf = basis_set_in%nsgf
     232        8820 :       ALLOCATE (basis_set_out%cgf_symbol(ncgf))
     233        8820 :       ALLOCATE (basis_set_out%sgf_symbol(nsgf))
     234       47662 :       basis_set_out%cgf_symbol = basis_set_in%cgf_symbol
     235       42508 :       basis_set_out%sgf_symbol = basis_set_in%sgf_symbol
     236        8820 :       ALLOCATE (basis_set_out%norm_cgf(ncgf))
     237       47662 :       basis_set_out%norm_cgf = basis_set_in%norm_cgf
     238        8820 :       ALLOCATE (basis_set_out%set_radius(nset))
     239       11826 :       basis_set_out%set_radius = basis_set_in%set_radius
     240       26460 :       ALLOCATE (basis_set_out%lmax(nset), basis_set_out%lmin(nset), basis_set_out%npgf(nset), basis_set_out%nshell(nset))
     241       11826 :       basis_set_out%lmax = basis_set_in%lmax
     242       11826 :       basis_set_out%lmin = basis_set_in%lmin
     243       11826 :       basis_set_out%npgf = basis_set_in%npgf
     244       11826 :       basis_set_out%nshell = basis_set_in%nshell
     245       26460 :       ALLOCATE (basis_set_out%lx(ncgf), basis_set_out%ly(ncgf), basis_set_out%lz(ncgf), basis_set_out%m(nsgf))
     246       47662 :       basis_set_out%lx = basis_set_in%lx
     247       47662 :       basis_set_out%ly = basis_set_in%ly
     248       47662 :       basis_set_out%lz = basis_set_in%lz
     249       42508 :       basis_set_out%m = basis_set_in%m
     250       14700 :       ALLOCATE (basis_set_out%ncgf_set(nset), basis_set_out%nsgf_set(nset))
     251       11826 :       basis_set_out%ncgf_set = basis_set_in%ncgf_set
     252       11826 :       basis_set_out%nsgf_set = basis_set_in%nsgf_set
     253        2940 :       maxco = SIZE(basis_set_in%cphi, 1)
     254       29400 :       ALLOCATE (basis_set_out%cphi(maxco, ncgf), basis_set_out%sphi(maxco, nsgf), basis_set_out%scon(maxco, nsgf))
     255     1246042 :       basis_set_out%cphi = basis_set_in%cphi
     256     1048940 :       basis_set_out%sphi = basis_set_in%sphi
     257     1048940 :       basis_set_out%scon = basis_set_in%scon
     258       11826 :       maxpgf = MAXVAL(basis_set_in%npgf)
     259       20580 :       ALLOCATE (basis_set_out%pgf_radius(maxpgf, nset), basis_set_out%zet(maxpgf, nset))
     260       44898 :       basis_set_out%pgf_radius = basis_set_in%pgf_radius
     261       44898 :       basis_set_out%zet = basis_set_in%zet
     262       11826 :       maxshell = MAXVAL(basis_set_in%nshell)
     263       20580 :       ALLOCATE (basis_set_out%first_cgf(maxshell, nset), basis_set_out%first_sgf(maxshell, nset))
     264       20580 :       ALLOCATE (basis_set_out%last_cgf(maxshell, nset), basis_set_out%last_sgf(maxshell, nset))
     265       32100 :       basis_set_out%first_cgf = basis_set_in%first_cgf
     266       32100 :       basis_set_out%first_sgf = basis_set_in%first_sgf
     267       32100 :       basis_set_out%last_cgf = basis_set_in%last_cgf
     268       32100 :       basis_set_out%last_sgf = basis_set_in%last_sgf
     269       20580 :       ALLOCATE (basis_set_out%n(maxshell, nset), basis_set_out%l(maxshell, nset))
     270       32100 :       basis_set_out%n = basis_set_in%n
     271       32100 :       basis_set_out%l = basis_set_in%l
     272       14700 :       ALLOCATE (basis_set_out%gcc(maxpgf, maxshell, nset))
     273      115370 :       basis_set_out%gcc = basis_set_in%gcc
     274             : 
     275        2940 :    END SUBROUTINE copy_gto_basis_set
     276             : 
     277             : ! **************************************************************************************************
     278             : !> \brief ...
     279             : !> \param basis_set ...
     280             : !> \param pbasis ...
     281             : ! **************************************************************************************************
     282          54 :    SUBROUTINE create_primitive_basis_set(basis_set, pbasis)
     283             : 
     284             :       ! Create a primitives only basis set
     285             : 
     286             :       TYPE(gto_basis_set_type), INTENT(IN)               :: basis_set
     287             :       TYPE(gto_basis_set_type), POINTER                  :: pbasis
     288             : 
     289             :       INTEGER                                            :: i, ico, ip, ipgf, iset, ishell, l, lm, &
     290             :                                                             lshell, m, maxco, mpgf, nc, ncgf, ns, &
     291             :                                                             nset, nsgf
     292          54 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: nindex, nprim
     293             :       REAL(KIND=dp)                                      :: zet0
     294          54 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: zet, zeta
     295             : 
     296         192 :       mpgf = SUM(basis_set%npgf)
     297         192 :       lm = MAXVAL(basis_set%lmax)
     298         594 :       ALLOCATE (zet(mpgf, 0:lm), zeta(mpgf, lm + 1), nindex(mpgf), nprim(0:lm))
     299        1784 :       zet = 0.0_dp
     300        1784 :       zeta = 0.0_dp
     301         196 :       DO l = 0, lm
     302         142 :          ip = 0
     303         560 :          DO iset = 1, basis_set%nset
     304         560 :             IF (basis_set%lmin(iset) <= l .AND. basis_set%lmax(iset) >= l) THEN
     305         634 :                DO ipgf = 1, basis_set%npgf(iset)
     306         492 :                   ip = ip + 1
     307         634 :                   zet(ip, l) = basis_set%zet(ipgf, iset)
     308             :                END DO
     309             :             END IF
     310             :          END DO
     311         196 :          nprim(l) = ip
     312             :       END DO
     313             : 
     314             :       ! sort exponents
     315         196 :       DO l = 0, lm
     316         634 :          zet(1:nprim(l), l) = -zet(1:nprim(l), l)
     317         142 :          CALL sort(zet(1:nprim(l), l), nprim(l), nindex)
     318             :          ! remove duplicates
     319         142 :          ip = 0
     320         142 :          zet0 = 0.0_dp
     321         634 :          DO i = 1, nprim(l)
     322         634 :             IF (ABS(zet0 - zet(i, l)) > 1.e-6_dp) THEN
     323         492 :                ip = ip + 1
     324         492 :                zeta(ip, l + 1) = zet(i, l)
     325             :             END IF
     326             :          END DO
     327         142 :          nprim(l) = ip
     328             :          !
     329         688 :          zeta(1:ip, l + 1) = -zeta(1:ip, l + 1)
     330             :       END DO
     331             : 
     332          54 :       CALL allocate_gto_basis_set(pbasis)
     333             : 
     334          54 :       pbasis%name = basis_set%name//"_primitive"
     335          54 :       pbasis%kind_radius = basis_set%kind_radius
     336          54 :       pbasis%short_kind_radius = basis_set%short_kind_radius
     337          54 :       pbasis%norm_type = basis_set%norm_type
     338          54 :       nset = lm + 1
     339          54 :       pbasis%nset = nset
     340         486 :       ALLOCATE (pbasis%lmax(nset), pbasis%lmin(nset), pbasis%npgf(nset), pbasis%nshell(nset))
     341         196 :       DO iset = 1, nset
     342         142 :          pbasis%lmax(iset) = iset - 1
     343         142 :          pbasis%lmin(iset) = iset - 1
     344         142 :          pbasis%npgf(iset) = nprim(iset - 1)
     345         196 :          pbasis%nshell(iset) = nprim(iset - 1)
     346             :       END DO
     347          54 :       pbasis%ncgf = 0
     348          54 :       pbasis%nsgf = 0
     349         196 :       DO l = 0, lm
     350         142 :          pbasis%ncgf = pbasis%ncgf + nprim(l)*((l + 1)*(l + 2))/2
     351         196 :          pbasis%nsgf = pbasis%nsgf + nprim(l)*(2*l + 1)
     352             :       END DO
     353         196 :       mpgf = MAXVAL(nprim)
     354         216 :       ALLOCATE (pbasis%zet(mpgf, nset))
     355         896 :       pbasis%zet(1:mpgf, 1:nset) = zeta(1:mpgf, 1:nset)
     356             : 
     357         378 :       ALLOCATE (pbasis%l(mpgf, nset), pbasis%n(mpgf, nset))
     358         196 :       DO iset = 1, nset
     359         688 :          DO ip = 1, nprim(iset - 1)
     360         492 :             pbasis%l(ip, iset) = iset - 1
     361         634 :             pbasis%n(ip, iset) = iset + ip - 1
     362             :          END DO
     363             :       END DO
     364             : 
     365         162 :       ALLOCATE (pbasis%cgf_symbol(pbasis%ncgf))
     366         162 :       ALLOCATE (pbasis%lx(pbasis%ncgf))
     367         162 :       ALLOCATE (pbasis%ly(pbasis%ncgf))
     368         162 :       ALLOCATE (pbasis%lz(pbasis%ncgf))
     369         162 :       ALLOCATE (pbasis%m(pbasis%nsgf))
     370         162 :       ALLOCATE (pbasis%sgf_symbol(pbasis%nsgf))
     371         270 :       ALLOCATE (pbasis%ncgf_set(nset), pbasis%nsgf_set(nset))
     372             : 
     373          54 :       ncgf = 0
     374          54 :       nsgf = 0
     375         196 :       DO iset = 1, nset
     376         142 :          l = iset - 1
     377         142 :          pbasis%ncgf_set(iset) = nprim(l)*((l + 1)*(l + 2))/2
     378         142 :          pbasis%nsgf_set(iset) = nprim(l)*(2*l + 1)
     379         688 :          DO ishell = 1, pbasis%nshell(iset)
     380         492 :             lshell = pbasis%l(ishell, iset)
     381        2458 :             DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
     382        1966 :                ncgf = ncgf + 1
     383        1966 :                pbasis%lx(ncgf) = indco(1, ico)
     384        1966 :                pbasis%ly(ncgf) = indco(2, ico)
     385        1966 :                pbasis%lz(ncgf) = indco(3, ico)
     386             :                pbasis%cgf_symbol(ncgf) = &
     387        8356 :                   cgf_symbol(pbasis%n(ishell, iset), (/pbasis%lx(ncgf), pbasis%ly(ncgf), pbasis%lz(ncgf)/))
     388             :             END DO
     389        2282 :             DO m = -lshell, lshell
     390        1648 :                nsgf = nsgf + 1
     391        1648 :                pbasis%m(nsgf) = m
     392        2140 :                pbasis%sgf_symbol(nsgf) = sgf_symbol(pbasis%n(ishell, iset), lshell, m)
     393             :             END DO
     394             :          END DO
     395             :       END DO
     396          54 :       CPASSERT(ncgf == pbasis%ncgf)
     397          54 :       CPASSERT(nsgf == pbasis%nsgf)
     398             : 
     399         270 :       ALLOCATE (pbasis%gcc(mpgf, mpgf, nset))
     400        5308 :       pbasis%gcc = 0.0_dp
     401         196 :       DO iset = 1, nset
     402         896 :          DO i = 1, mpgf
     403         842 :             pbasis%gcc(i, i, iset) = 1.0_dp
     404             :          END DO
     405             :       END DO
     406             : 
     407         216 :       ALLOCATE (pbasis%first_cgf(mpgf, nset))
     408         216 :       ALLOCATE (pbasis%first_sgf(mpgf, nset))
     409         216 :       ALLOCATE (pbasis%last_cgf(mpgf, nset))
     410         216 :       ALLOCATE (pbasis%last_sgf(mpgf, nset))
     411          54 :       nc = 0
     412          54 :       ns = 0
     413          54 :       maxco = 0
     414         196 :       DO iset = 1, nset
     415         634 :          DO ishell = 1, pbasis%nshell(iset)
     416         492 :             lshell = pbasis%l(ishell, iset)
     417         492 :             pbasis%first_cgf(ishell, iset) = nc + 1
     418         492 :             nc = nc + nco(lshell)
     419         492 :             pbasis%last_cgf(ishell, iset) = nc
     420         492 :             pbasis%first_sgf(ishell, iset) = ns + 1
     421         492 :             ns = ns + nso(lshell)
     422         634 :             pbasis%last_sgf(ishell, iset) = ns
     423             :          END DO
     424         196 :          maxco = MAX(maxco, pbasis%npgf(iset)*ncoset(pbasis%lmax(iset)))
     425             :       END DO
     426             : 
     427         162 :       ALLOCATE (pbasis%norm_cgf(ncgf))
     428         216 :       ALLOCATE (pbasis%cphi(maxco, ncgf))
     429      171472 :       pbasis%cphi = 0.0_dp
     430         216 :       ALLOCATE (pbasis%sphi(maxco, nsgf))
     431      136270 :       pbasis%sphi = 0.0_dp
     432         216 :       ALLOCATE (pbasis%scon(maxco, ncgf))
     433      171472 :       pbasis%scon = 0.0_dp
     434         162 :       ALLOCATE (pbasis%set_radius(nset))
     435         216 :       ALLOCATE (pbasis%pgf_radius(mpgf, nset))
     436         896 :       pbasis%pgf_radius = 0.0_dp
     437             : 
     438          54 :       CALL init_orb_basis_set(pbasis)
     439             : 
     440          54 :       DEALLOCATE (zet, zeta, nindex, nprim)
     441             : 
     442          54 :    END SUBROUTINE create_primitive_basis_set
     443             : 
     444             : ! **************************************************************************************************
     445             : !> \brief ...
     446             : !> \param basis_set ...
     447             : !> \param basis_set_add ...
     448             : ! **************************************************************************************************
     449          42 :    SUBROUTINE combine_basis_sets(basis_set, basis_set_add)
     450             : 
     451             :       ! Combine two Gaussian-type orbital (GTO) basis sets.
     452             : 
     453             :       TYPE(gto_basis_set_type), INTENT(INOUT)            :: basis_set
     454             :       TYPE(gto_basis_set_type), INTENT(IN)               :: basis_set_add
     455             : 
     456          42 :       CHARACTER(LEN=12), DIMENSION(:), POINTER           :: cgf_symbol
     457          42 :       CHARACTER(LEN=6), DIMENSION(:), POINTER            :: sgf_symbol
     458             :       INTEGER                                            :: iset, ishell, lshell, maxco, maxpgf, &
     459             :                                                             maxshell, nc, ncgf, ncgfn, ncgfo, ns, &
     460             :                                                             nset, nsetn, nseto, nsgf, nsgfn, nsgfo
     461             : 
     462          84 :       basis_set%name = basis_set%name//basis_set_add%name
     463          42 :       basis_set%nset = basis_set%nset + basis_set_add%nset
     464          42 :       basis_set%ncgf = basis_set%ncgf + basis_set_add%ncgf
     465          42 :       basis_set%nsgf = basis_set%nsgf + basis_set_add%nsgf
     466          42 :       nset = basis_set%nset
     467          42 :       ncgf = basis_set%ncgf
     468          42 :       nsgf = basis_set%nsgf
     469             : 
     470          42 :       nsetn = basis_set_add%nset
     471          42 :       nseto = nset - nsetn
     472          42 :       CALL reallocate(basis_set%set_radius, 1, nset) ! to be defined later
     473          42 :       CALL reallocate(basis_set%lmax, 1, nset)
     474          42 :       CALL reallocate(basis_set%lmin, 1, nset)
     475          42 :       CALL reallocate(basis_set%npgf, 1, nset)
     476          42 :       CALL reallocate(basis_set%nshell, 1, nset)
     477         160 :       basis_set%lmax(nseto + 1:nset) = basis_set_add%lmax(1:nsetn)
     478         160 :       basis_set%lmin(nseto + 1:nset) = basis_set_add%lmin(1:nsetn)
     479         160 :       basis_set%npgf(nseto + 1:nset) = basis_set_add%npgf(1:nsetn)
     480         160 :       basis_set%nshell(nseto + 1:nset) = basis_set_add%nshell(1:nsetn)
     481          42 :       CALL reallocate(basis_set%ncgf_set, 1, nset)
     482          42 :       CALL reallocate(basis_set%nsgf_set, 1, nset)
     483         160 :       basis_set%ncgf_set(nseto + 1:nset) = basis_set_add%ncgf_set(1:nsetn)
     484         160 :       basis_set%nsgf_set(nseto + 1:nset) = basis_set_add%nsgf_set(1:nsetn)
     485             : 
     486          42 :       nsgfn = basis_set_add%nsgf
     487          42 :       nsgfo = nsgf - nsgfn
     488          42 :       ncgfn = basis_set_add%ncgf
     489          42 :       ncgfo = ncgf - ncgfn
     490             : 
     491         210 :       ALLOCATE (cgf_symbol(ncgf), sgf_symbol(nsgf))
     492         562 :       cgf_symbol(1:ncgfo) = basis_set%cgf_symbol(1:ncgfo)
     493        1868 :       cgf_symbol(ncgfo + 1:ncgf) = basis_set_add%cgf_symbol(1:ncgfn)
     494         530 :       sgf_symbol(1:nsgfo) = basis_set%sgf_symbol(1:nsgfo)
     495        1554 :       sgf_symbol(nsgfo + 1:nsgf) = basis_set_add%sgf_symbol(1:nsgfn)
     496          42 :       DEALLOCATE (basis_set%cgf_symbol, basis_set%sgf_symbol)
     497         210 :       ALLOCATE (basis_set%cgf_symbol(ncgf), basis_set%sgf_symbol(nsgf))
     498        2388 :       basis_set%cgf_symbol = cgf_symbol
     499        2042 :       basis_set%sgf_symbol = sgf_symbol
     500          42 :       DEALLOCATE (cgf_symbol, sgf_symbol)
     501             : 
     502          42 :       CALL reallocate(basis_set%lx, 1, ncgf)
     503          42 :       CALL reallocate(basis_set%ly, 1, ncgf)
     504          42 :       CALL reallocate(basis_set%lz, 1, ncgf)
     505          42 :       CALL reallocate(basis_set%m, 1, nsgf)
     506        1868 :       basis_set%lx(ncgfo + 1:ncgf) = basis_set_add%lx(1:ncgfn)
     507        1868 :       basis_set%ly(ncgfo + 1:ncgf) = basis_set_add%ly(1:ncgfn)
     508        1868 :       basis_set%lz(ncgfo + 1:ncgf) = basis_set_add%lz(1:ncgfn)
     509        1554 :       basis_set%m(nsgfo + 1:nsgf) = basis_set_add%m(1:nsgfn)
     510             : 
     511         238 :       maxpgf = MAXVAL(basis_set%npgf)
     512          42 :       CALL reallocate(basis_set%zet, 1, maxpgf, 1, nset)
     513          42 :       nc = SIZE(basis_set_add%zet, 1)
     514         160 :       DO iset = 1, nsetn
     515         770 :          basis_set%zet(1:nc, nseto + iset) = basis_set_add%zet(1:nc, iset)
     516             :       END DO
     517             : 
     518         238 :       maxshell = MAXVAL(basis_set%nshell)
     519          42 :       CALL reallocate(basis_set%l, 1, maxshell, 1, nset)
     520          42 :       CALL reallocate(basis_set%n, 1, maxshell, 1, nset)
     521          42 :       nc = SIZE(basis_set_add%l, 1)
     522         160 :       DO iset = 1, nsetn
     523         728 :          basis_set%l(1:nc, nseto + iset) = basis_set_add%l(1:nc, iset)
     524         770 :          basis_set%n(1:nc, nseto + iset) = basis_set_add%n(1:nc, iset)
     525             :       END DO
     526             : 
     527          42 :       CALL reallocate(basis_set%first_cgf, 1, maxshell, 1, nset)
     528          42 :       CALL reallocate(basis_set%first_sgf, 1, maxshell, 1, nset)
     529          42 :       CALL reallocate(basis_set%last_cgf, 1, maxshell, 1, nset)
     530          42 :       CALL reallocate(basis_set%last_sgf, 1, maxshell, 1, nset)
     531          42 :       nc = 0
     532          42 :       ns = 0
     533         238 :       DO iset = 1, nset
     534         866 :          DO ishell = 1, basis_set%nshell(iset)
     535         628 :             lshell = basis_set%l(ishell, iset)
     536         628 :             basis_set%first_cgf(ishell, iset) = nc + 1
     537         628 :             nc = nc + nco(lshell)
     538         628 :             basis_set%last_cgf(ishell, iset) = nc
     539         628 :             basis_set%first_sgf(ishell, iset) = ns + 1
     540         628 :             ns = ns + nso(lshell)
     541         824 :             basis_set%last_sgf(ishell, iset) = ns
     542             :          END DO
     543             :       END DO
     544             : 
     545          42 :       CALL reallocate(basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)
     546          42 :       nc = SIZE(basis_set_add%gcc, 1)
     547          42 :       ns = SIZE(basis_set_add%gcc, 2)
     548         160 :       DO iset = 1, nsetn
     549        4840 :          basis_set%gcc(1:nc, 1:ns, nseto + iset) = basis_set_add%gcc(1:nc, 1:ns, iset)
     550             :       END DO
     551             : 
     552             :       ! these arrays are determined later using initialization calls
     553          42 :       CALL reallocate(basis_set%norm_cgf, 1, ncgf)
     554          42 :       maxco = MAX(SIZE(basis_set%cphi, 1), SIZE(basis_set_add%cphi, 1))
     555          42 :       CALL reallocate(basis_set%cphi, 1, maxco, 1, ncgf)
     556          42 :       CALL reallocate(basis_set%sphi, 1, maxco, 1, nsgf)
     557          42 :       CALL reallocate(basis_set%scon, 1, maxco, 1, nsgf)
     558          42 :       CALL reallocate(basis_set%pgf_radius, 1, maxpgf, 1, nset)
     559             : 
     560          42 :    END SUBROUTINE combine_basis_sets
     561             : 
     562             : ! **************************************************************************************************
     563             : !> \brief ...
     564             : !> \param gto_basis_set ...
     565             : !> \param name ...
     566             : !> \param aliases ...
     567             : !> \param norm_type ...
     568             : !> \param kind_radius ...
     569             : !> \param ncgf ...
     570             : !> \param nset ...
     571             : !> \param nsgf ...
     572             : !> \param cgf_symbol ...
     573             : !> \param sgf_symbol ...
     574             : !> \param norm_cgf ...
     575             : !> \param set_radius ...
     576             : !> \param lmax ...
     577             : !> \param lmin ...
     578             : !> \param lx ...
     579             : !> \param ly ...
     580             : !> \param lz ...
     581             : !> \param m ...
     582             : !> \param ncgf_set ...
     583             : !> \param npgf ...
     584             : !> \param nsgf_set ...
     585             : !> \param nshell ...
     586             : !> \param cphi ...
     587             : !> \param pgf_radius ...
     588             : !> \param sphi ...
     589             : !> \param scon ...
     590             : !> \param zet ...
     591             : !> \param first_cgf ...
     592             : !> \param first_sgf ...
     593             : !> \param l ...
     594             : !> \param last_cgf ...
     595             : !> \param last_sgf ...
     596             : !> \param n ...
     597             : !> \param gcc ...
     598             : !> \param maxco ...
     599             : !> \param maxl ...
     600             : !> \param maxpgf ...
     601             : !> \param maxsgf_set ...
     602             : !> \param maxshell ...
     603             : !> \param maxso ...
     604             : !> \param nco_sum ...
     605             : !> \param npgf_sum ...
     606             : !> \param nshell_sum ...
     607             : !> \param maxder ...
     608             : !> \param short_kind_radius ...
     609             : ! **************************************************************************************************
     610    32115071 :    SUBROUTINE get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, &
     611             :                                 nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, &
     612             :                                 m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, &
     613             :                                 last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, &
     614             :                                 npgf_sum, nshell_sum, maxder, short_kind_radius)
     615             : 
     616             :       ! Get informations about a Gaussian-type orbital (GTO) basis set.
     617             : 
     618             :       ! - Creation (10.01.2002,MK)
     619             : 
     620             :       TYPE(gto_basis_set_type), INTENT(IN)               :: gto_basis_set
     621             :       CHARACTER(LEN=default_string_length), &
     622             :          INTENT(OUT), OPTIONAL                           :: name, aliases
     623             :       INTEGER, INTENT(OUT), OPTIONAL                     :: norm_type
     624             :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: kind_radius
     625             :       INTEGER, INTENT(OUT), OPTIONAL                     :: ncgf, nset, nsgf
     626             :       CHARACTER(LEN=12), DIMENSION(:), OPTIONAL, POINTER :: cgf_symbol
     627             :       CHARACTER(LEN=6), DIMENSION(:), OPTIONAL, POINTER  :: sgf_symbol
     628             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: norm_cgf, set_radius
     629             :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: lmax, lmin, lx, ly, lz, m, ncgf_set, &
     630             :                                                             npgf, nsgf_set, nshell
     631             :       REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: cphi, pgf_radius, sphi, scon, zet
     632             :       INTEGER, DIMENSION(:, :), OPTIONAL, POINTER        :: first_cgf, first_sgf, l, last_cgf, &
     633             :                                                             last_sgf, n
     634             :       REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
     635             :          POINTER                                         :: gcc
     636             :       INTEGER, INTENT(OUT), OPTIONAL                     :: maxco, maxl, maxpgf, maxsgf_set, &
     637             :                                                             maxshell, maxso, nco_sum, npgf_sum, &
     638             :                                                             nshell_sum
     639             :       INTEGER, INTENT(IN), OPTIONAL                      :: maxder
     640             :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: short_kind_radius
     641             : 
     642             :       INTEGER                                            :: iset, nder
     643             : 
     644    32115071 :       IF (PRESENT(name)) name = gto_basis_set%name
     645    32115071 :       IF (PRESENT(aliases)) aliases = gto_basis_set%aliases
     646    32115071 :       IF (PRESENT(norm_type)) norm_type = gto_basis_set%norm_type
     647    32115071 :       IF (PRESENT(kind_radius)) kind_radius = gto_basis_set%kind_radius
     648    32115071 :       IF (PRESENT(short_kind_radius)) short_kind_radius = gto_basis_set%short_kind_radius
     649    32115071 :       IF (PRESENT(ncgf)) ncgf = gto_basis_set%ncgf
     650    32115071 :       IF (PRESENT(nset)) nset = gto_basis_set%nset
     651    32115071 :       IF (PRESENT(nsgf)) nsgf = gto_basis_set%nsgf
     652    32115071 :       IF (PRESENT(cgf_symbol)) cgf_symbol => gto_basis_set%cgf_symbol
     653    32115071 :       IF (PRESENT(sgf_symbol)) sgf_symbol => gto_basis_set%sgf_symbol
     654    32115071 :       IF (PRESENT(norm_cgf)) norm_cgf => gto_basis_set%norm_cgf
     655    32115071 :       IF (PRESENT(set_radius)) set_radius => gto_basis_set%set_radius
     656    32115071 :       IF (PRESENT(lmax)) lmax => gto_basis_set%lmax
     657    32115071 :       IF (PRESENT(lmin)) lmin => gto_basis_set%lmin
     658    32115071 :       IF (PRESENT(lx)) lx => gto_basis_set%lx
     659    32115071 :       IF (PRESENT(ly)) ly => gto_basis_set%ly
     660    32115071 :       IF (PRESENT(lz)) lz => gto_basis_set%lz
     661    32115071 :       IF (PRESENT(m)) m => gto_basis_set%m
     662    32115071 :       IF (PRESENT(ncgf_set)) ncgf_set => gto_basis_set%ncgf_set
     663    32115071 :       IF (PRESENT(npgf)) npgf => gto_basis_set%npgf
     664    32115071 :       IF (PRESENT(nsgf_set)) nsgf_set => gto_basis_set%nsgf_set
     665    32115071 :       IF (PRESENT(nshell)) nshell => gto_basis_set%nshell
     666    32115071 :       IF (PRESENT(cphi)) cphi => gto_basis_set%cphi
     667    32115071 :       IF (PRESENT(pgf_radius)) pgf_radius => gto_basis_set%pgf_radius
     668    32115071 :       IF (PRESENT(sphi)) sphi => gto_basis_set%sphi
     669    32115071 :       IF (PRESENT(scon)) scon => gto_basis_set%scon
     670    32115071 :       IF (PRESENT(zet)) zet => gto_basis_set%zet
     671    32115071 :       IF (PRESENT(first_cgf)) first_cgf => gto_basis_set%first_cgf
     672    32115071 :       IF (PRESENT(first_sgf)) first_sgf => gto_basis_set%first_sgf
     673    32115071 :       IF (PRESENT(l)) l => gto_basis_set%l
     674    32115071 :       IF (PRESENT(last_cgf)) last_cgf => gto_basis_set%last_cgf
     675    32115071 :       IF (PRESENT(last_sgf)) last_sgf => gto_basis_set%last_sgf
     676    32115071 :       IF (PRESENT(n)) n => gto_basis_set%n
     677    32115071 :       IF (PRESENT(gcc)) gcc => gto_basis_set%gcc
     678    32115071 :       IF (PRESENT(maxco)) THEN
     679     5336753 :          maxco = 0
     680     5336753 :          IF (PRESENT(maxder)) THEN
     681           0 :             nder = maxder
     682             :          ELSE
     683             :             nder = 0
     684             :          END IF
     685    17169816 :          DO iset = 1, gto_basis_set%nset
     686             :             maxco = MAX(maxco, gto_basis_set%npgf(iset)* &
     687    17169816 :                         ncoset(gto_basis_set%lmax(iset) + nder))
     688             :          END DO
     689             :       END IF
     690    32115071 :       IF (PRESENT(maxl)) THEN
     691     5084041 :          maxl = -1
     692    15611364 :          DO iset = 1, gto_basis_set%nset
     693    15611364 :             maxl = MAX(maxl, gto_basis_set%lmax(iset))
     694             :          END DO
     695             :       END IF
     696    32115071 :       IF (PRESENT(maxpgf)) THEN
     697        4222 :          maxpgf = 0
     698       18662 :          DO iset = 1, gto_basis_set%nset
     699       18662 :             maxpgf = MAX(maxpgf, gto_basis_set%npgf(iset))
     700             :          END DO
     701             :       END IF
     702    32115071 :       IF (PRESENT(maxsgf_set)) THEN
     703      428861 :          maxsgf_set = 0
     704     2161540 :          DO iset = 1, gto_basis_set%nset
     705     2161540 :             maxsgf_set = MAX(maxsgf_set, gto_basis_set%nsgf_set(iset))
     706             :          END DO
     707             :       END IF
     708    32115071 :       IF (PRESENT(maxshell)) THEN ! MAXVAL on structure component avoided
     709        2608 :          maxshell = 0
     710       11902 :          DO iset = 1, gto_basis_set%nset
     711       11902 :             maxshell = MAX(maxshell, gto_basis_set%nshell(iset))
     712             :          END DO
     713             :       END IF
     714    32115071 :       IF (PRESENT(maxso)) THEN
     715     1791928 :          maxso = 0
     716     6580900 :          DO iset = 1, gto_basis_set%nset
     717             :             maxso = MAX(maxso, gto_basis_set%npgf(iset)* &
     718     6580900 :                         nsoset(gto_basis_set%lmax(iset)))
     719             :          END DO
     720             :       END IF
     721             : 
     722    32115071 :       IF (PRESENT(nco_sum)) THEN
     723       73036 :          nco_sum = 0
     724      257854 :          DO iset = 1, gto_basis_set%nset
     725             :             nco_sum = nco_sum + gto_basis_set%npgf(iset)* &
     726      257854 :                       ncoset(gto_basis_set%lmax(iset))
     727             :          END DO
     728             :       END IF
     729    32128993 :       IF (PRESENT(npgf_sum)) npgf_sum = SUM(gto_basis_set%npgf)
     730    32128993 :       IF (PRESENT(nshell_sum)) nshell_sum = SUM(gto_basis_set%nshell)
     731             : 
     732    32115071 :    END SUBROUTINE get_gto_basis_set
     733             : 
     734             : ! **************************************************************************************************
     735             : !> \brief ...
     736             : !> \param gto_basis_set ...
     737             : ! **************************************************************************************************
     738          12 :    SUBROUTINE init_aux_basis_set(gto_basis_set)
     739             : 
     740             :       ! Initialise a Gaussian-type orbital (GTO) basis set data set.
     741             : 
     742             :       ! - Creation (06.12.2000,MK)
     743             : 
     744             :       TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set
     745             : 
     746             :       CHARACTER(len=*), PARAMETER :: routineN = 'init_aux_basis_set'
     747             : 
     748             :       INTEGER                                            :: handle
     749             : 
     750             : ! -------------------------------------------------------------------------
     751             : 
     752           6 :       IF (.NOT. ASSOCIATED(gto_basis_set)) RETURN
     753             : 
     754           6 :       CALL timeset(routineN, handle)
     755             : 
     756          12 :       SELECT CASE (gto_basis_set%norm_type)
     757             :       CASE (0)
     758             :          ! No normalisation requested
     759             :       CASE (1)
     760           6 :          CALL init_norm_cgf_aux_2(gto_basis_set)
     761             :       CASE (2)
     762             :          ! WARNING this was never tested
     763           0 :          CALL init_norm_cgf_aux(gto_basis_set)
     764             :       CASE DEFAULT
     765           6 :          CPABORT("Normalization method not specified")
     766             :       END SELECT
     767             : 
     768             :       ! Initialise the transformation matrices "pgf" -> "cgf"
     769           6 :       CALL init_cphi_and_sphi(gto_basis_set)
     770             : 
     771           6 :       CALL timestop(handle)
     772             : 
     773             :    END SUBROUTINE init_aux_basis_set
     774             : 
     775             : ! **************************************************************************************************
     776             : !> \brief ...
     777             : !> \param gto_basis_set ...
     778             : ! **************************************************************************************************
     779       16531 :    SUBROUTINE init_cphi_and_sphi(gto_basis_set)
     780             : 
     781             :       ! Initialise the matrices for the transformation of primitive Cartesian
     782             :       ! Gaussian-type functions to contracted Cartesian (cphi) and spherical
     783             :       ! (sphi) Gaussian-type functions.
     784             : 
     785             :       ! - Creation (20.09.2000,MK)
     786             : 
     787             :       TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set
     788             : 
     789             :       CHARACTER(len=*), PARAMETER :: routineN = 'init_cphi_and_sphi'
     790             : 
     791             :       INTEGER                                            :: first_cgf, first_sgf, handle, icgf, ico, &
     792             :                                                             ipgf, iset, ishell, l, last_sgf, lmax, &
     793             :                                                             lmin, n, n1, n2, ncgf, nn, nn1, nn2, &
     794             :                                                             npgf, nsgf
     795             : 
     796             : ! -------------------------------------------------------------------------
     797             : ! Build the Cartesian transformation matrix "cphi"
     798             : 
     799       16531 :       CALL timeset(routineN, handle)
     800             : 
     801     6354740 :       gto_basis_set%cphi = 0.0_dp
     802       61272 :       DO iset = 1, gto_basis_set%nset
     803       44741 :          n = ncoset(gto_basis_set%lmax(iset))
     804      135201 :          DO ishell = 1, gto_basis_set%nshell(iset)
     805      274265 :             DO icgf = gto_basis_set%first_cgf(ishell, iset), &
     806      118670 :                gto_basis_set%last_cgf(ishell, iset)
     807             :                ico = coset(gto_basis_set%lx(icgf), &
     808             :                            gto_basis_set%ly(icgf), &
     809      200336 :                            gto_basis_set%lz(icgf))
     810      901577 :                DO ipgf = 1, gto_basis_set%npgf(iset)
     811             :                   gto_basis_set%cphi(ico, icgf) = gto_basis_set%norm_cgf(icgf)* &
     812      627312 :                                                   gto_basis_set%gcc(ipgf, ishell, iset)
     813      827648 :                   ico = ico + n
     814             :                END DO
     815             :             END DO
     816             :          END DO
     817             :       END DO
     818             : 
     819             :       ! Build the spherical transformation matrix "sphi"
     820             : 
     821       16531 :       n = SIZE(gto_basis_set%cphi, 1)
     822             : 
     823     5552713 :       gto_basis_set%sphi = 0.0_dp
     824       16531 :       IF (n .GT. 0) THEN
     825       61256 :          DO iset = 1, gto_basis_set%nset
     826      135179 :             DO ishell = 1, gto_basis_set%nshell(iset)
     827       73923 :                l = gto_basis_set%l(ishell, iset)
     828       73923 :                first_cgf = gto_basis_set%first_cgf(ishell, iset)
     829       73923 :                first_sgf = gto_basis_set%first_sgf(ishell, iset)
     830       73923 :                ncgf = nco(l)
     831       73923 :                nsgf = nso(l)
     832             :                CALL dgemm("N", "T", n, nsgf, ncgf, &
     833             :                           1.0_dp, gto_basis_set%cphi(1, first_cgf), n, &
     834             :                           orbtramat(l)%c2s(1, 1), nsgf, &
     835      118658 :                           0.0_dp, gto_basis_set%sphi(1, first_sgf), n)
     836             :             END DO
     837             :          END DO
     838             :       END IF
     839             : 
     840             :       ! Build the reduced transformation matrix "scon"
     841             :       ! This matrix transforms from cartesian primitifs to spherical contracted functions
     842             :       ! "scon" only includes primitifs (lmin -> lmax), whereas sphi is (0 -> lmax)
     843             : 
     844       16531 :       n = SIZE(gto_basis_set%scon, 1)
     845             : 
     846     5587915 :       gto_basis_set%scon = 0.0_dp
     847       16531 :       IF (n .GT. 0) THEN
     848       61256 :          DO iset = 1, gto_basis_set%nset
     849       44735 :             lmin = gto_basis_set%lmin(iset)
     850       44735 :             lmax = gto_basis_set%lmax(iset)
     851       44735 :             npgf = gto_basis_set%npgf(iset)
     852       44735 :             nn = ncoset(lmax) - ncoset(lmin - 1)
     853      135179 :             DO ishell = 1, gto_basis_set%nshell(iset)
     854       73923 :                first_sgf = gto_basis_set%first_sgf(ishell, iset)
     855       73923 :                last_sgf = gto_basis_set%last_sgf(ishell, iset)
     856      381779 :                DO ipgf = 1, npgf
     857      263121 :                   nn1 = (ipgf - 1)*ncoset(lmax) + ncoset(lmin - 1) + 1
     858      263121 :                   nn2 = ipgf*ncoset(lmax)
     859      263121 :                   n1 = (ipgf - 1)*nn + 1
     860      263121 :                   n2 = ipgf*nn
     861     4711621 :                   gto_basis_set%scon(n1:n2, first_sgf:last_sgf) = gto_basis_set%sphi(nn1:nn2, first_sgf:last_sgf)
     862             :                END DO
     863             :             END DO
     864             :          END DO
     865             :       END IF
     866             : 
     867       16531 :       CALL timestop(handle)
     868             : 
     869       16531 :    END SUBROUTINE init_cphi_and_sphi
     870             : 
     871             : ! **************************************************************************************************
     872             : !> \brief ...
     873             : !> \param gto_basis_set ...
     874             : ! **************************************************************************************************
     875           0 :    SUBROUTINE init_norm_cgf_aux(gto_basis_set)
     876             : 
     877             :       ! Initialise the normalization factors of the contracted Cartesian Gaussian
     878             :       ! functions, if the Gaussian functions represent charge distributions.
     879             : 
     880             :       ! - Creation (07.12.2000,MK)
     881             : 
     882             :       TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set
     883             : 
     884             :       INTEGER                                            :: icgf, ico, ipgf, iset, ishell, jco, &
     885             :                                                             jpgf, ll, lmax, lmin, lx, ly, lz, n, &
     886             :                                                             npgfa
     887             :       REAL(KIND=dp)                                      :: fnorm, gcca, gccb
     888           0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: ff
     889             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: gaa
     890             :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :)     :: vv
     891           0 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: rpgfa, zeta
     892             : 
     893             : ! -------------------------------------------------------------------------
     894             : 
     895           0 :       n = 0
     896           0 :       ll = 0
     897           0 :       DO iset = 1, gto_basis_set%nset
     898           0 :          n = MAX(n, gto_basis_set%npgf(iset)*ncoset(gto_basis_set%lmax(iset)))
     899           0 :          ll = MAX(ll, gto_basis_set%lmax(iset))
     900             :       END DO
     901             : 
     902           0 :       ALLOCATE (gaa(n, n))
     903           0 :       ALLOCATE (vv(ncoset(ll), ncoset(ll), ll + ll + 1))
     904           0 :       ALLOCATE (ff(0:ll + ll))
     905             : 
     906           0 :       DO iset = 1, gto_basis_set%nset
     907           0 :          lmax = gto_basis_set%lmax(iset)
     908           0 :          lmin = gto_basis_set%lmin(iset)
     909           0 :          n = ncoset(lmax)
     910           0 :          npgfa = gto_basis_set%npgf(iset)
     911           0 :          rpgfa => gto_basis_set%pgf_radius(1:npgfa, iset)
     912           0 :          zeta => gto_basis_set%zet(1:npgfa, iset)
     913             :          CALL coulomb2(lmax, npgfa, zeta, rpgfa, lmin, &
     914             :                        lmax, npgfa, zeta, rpgfa, lmin, &
     915           0 :                        (/0.0_dp, 0.0_dp, 0.0_dp/), 0.0_dp, gaa, vv, ff(0:))
     916           0 :          DO ishell = 1, gto_basis_set%nshell(iset)
     917           0 :             DO icgf = gto_basis_set%first_cgf(ishell, iset), &
     918           0 :                gto_basis_set%last_cgf(ishell, iset)
     919           0 :                lx = gto_basis_set%lx(icgf)
     920           0 :                ly = gto_basis_set%ly(icgf)
     921           0 :                lz = gto_basis_set%lz(icgf)
     922           0 :                ico = coset(lx, ly, lz)
     923           0 :                fnorm = 0.0_dp
     924           0 :                DO ipgf = 1, npgfa
     925           0 :                   gcca = gto_basis_set%gcc(ipgf, ishell, iset)
     926           0 :                   jco = coset(lx, ly, lz)
     927           0 :                   DO jpgf = 1, npgfa
     928           0 :                      gccb = gto_basis_set%gcc(jpgf, ishell, iset)
     929           0 :                      fnorm = fnorm + gcca*gccb*gaa(ico, jco)
     930           0 :                      jco = jco + n
     931             :                   END DO
     932           0 :                   ico = ico + n
     933             :                END DO
     934           0 :                gto_basis_set%norm_cgf(icgf) = 1.0_dp/SQRT(fnorm)
     935             :             END DO
     936             :          END DO
     937             :       END DO
     938             : 
     939           0 :       DEALLOCATE (vv, ff)
     940           0 :       DEALLOCATE (gaa)
     941             : 
     942           0 :    END SUBROUTINE init_norm_cgf_aux
     943             : 
     944             : ! **************************************************************************************************
     945             : !> \brief ...
     946             : !> \param gto_basis_set ...
     947             : ! **************************************************************************************************
     948           6 :    ELEMENTAL SUBROUTINE init_norm_cgf_aux_2(gto_basis_set)
     949             : 
     950             :       ! Initialise the normalization factors of the auxiliary Cartesian Gaussian
     951             :       ! functions (Kim-Gordon polarization basis) Norm = 1.
     952             : 
     953             :       ! - Creation (07.12.2000,GT)
     954             : 
     955             :       TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set
     956             : 
     957             :       INTEGER                                            :: icgf, iset, ishell
     958             : 
     959          92 :       DO iset = 1, gto_basis_set%nset
     960         178 :          DO ishell = 1, gto_basis_set%nshell(iset)
     961         396 :             DO icgf = gto_basis_set%first_cgf(ishell, iset), &
     962         172 :                gto_basis_set%last_cgf(ishell, iset)
     963         396 :                gto_basis_set%norm_cgf(icgf) = 1.0_dp
     964             :             END DO
     965             :          END DO
     966             :       END DO
     967             : 
     968           6 :    END SUBROUTINE init_norm_cgf_aux_2
     969             : 
     970             : ! **************************************************************************************************
     971             : !> \brief Initialise the normalization factors of the contracted Cartesian Gaussian functions.
     972             : !> \param gto_basis_set ...
     973             : !> \author MK
     974             : ! **************************************************************************************************
     975       14911 :    ELEMENTAL SUBROUTINE init_norm_cgf_orb(gto_basis_set)
     976             : 
     977             :       TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set
     978             : 
     979             :       INTEGER                                            :: icgf, ipgf, iset, ishell, jpgf, l, lx, &
     980             :                                                             ly, lz
     981             :       REAL(KIND=dp)                                      :: expzet, fnorm, gcca, gccb, prefac, zeta, &
     982             :                                                             zetb
     983             : 
     984       54496 :       DO iset = 1, gto_basis_set%nset
     985      120325 :          DO ishell = 1, gto_basis_set%nshell(iset)
     986             : 
     987       65829 :             l = gto_basis_set%l(ishell, iset)
     988             : 
     989       65829 :             expzet = 0.5_dp*REAL(2*l + 3, dp)
     990             : 
     991       65829 :             fnorm = 0.0_dp
     992             : 
     993      313334 :             DO ipgf = 1, gto_basis_set%npgf(iset)
     994      247505 :                gcca = gto_basis_set%gcc(ipgf, ishell, iset)
     995      247505 :                zeta = gto_basis_set%zet(ipgf, iset)
     996     1801745 :                DO jpgf = 1, gto_basis_set%npgf(iset)
     997     1488411 :                   gccb = gto_basis_set%gcc(jpgf, ishell, iset)
     998     1488411 :                   zetb = gto_basis_set%zet(jpgf, iset)
     999     1735916 :                   fnorm = fnorm + gcca*gccb/(zeta + zetb)**expzet
    1000             :                END DO
    1001             :             END DO
    1002             : 
    1003       65829 :             fnorm = 0.5_dp**l*pi**1.5_dp*fnorm
    1004             : 
    1005      246887 :             DO icgf = gto_basis_set%first_cgf(ishell, iset), &
    1006      105414 :                gto_basis_set%last_cgf(ishell, iset)
    1007      181058 :                lx = gto_basis_set%lx(icgf)
    1008      181058 :                ly = gto_basis_set%ly(icgf)
    1009      181058 :                lz = gto_basis_set%lz(icgf)
    1010      181058 :                prefac = dfac(2*lx - 1)*dfac(2*ly - 1)*dfac(2*lz - 1)
    1011      246887 :                gto_basis_set%norm_cgf(icgf) = 1.0_dp/SQRT(prefac*fnorm)
    1012             :             END DO
    1013             : 
    1014             :          END DO
    1015             :       END DO
    1016             : 
    1017       14911 :    END SUBROUTINE init_norm_cgf_orb
    1018             : 
    1019             : ! **************************************************************************************************
    1020             : !> \brief Initialise the normalization factors of the contracted Cartesian Gaussian
    1021             : !>        functions used for frozen density representation.
    1022             : !> \param gto_basis_set ...
    1023             : !> \author GT
    1024             : ! **************************************************************************************************
    1025           0 :    ELEMENTAL SUBROUTINE init_norm_cgf_orb_den(gto_basis_set)
    1026             : 
    1027             :       TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set
    1028             : 
    1029             :       INTEGER                                            :: icgf, ipgf, iset, ishell, l
    1030             :       REAL(KIND=dp)                                      :: expzet, gcca, prefac, zeta
    1031             : 
    1032           0 :       DO iset = 1, gto_basis_set%nset
    1033           0 :          DO ishell = 1, gto_basis_set%nshell(iset)
    1034           0 :             l = gto_basis_set%l(ishell, iset)
    1035           0 :             expzet = 0.5_dp*REAL(2*l + 3, dp)
    1036           0 :             prefac = (1.0_dp/pi)**1.5_dp
    1037           0 :             DO ipgf = 1, gto_basis_set%npgf(iset)
    1038           0 :                gcca = gto_basis_set%gcc(ipgf, ishell, iset)
    1039           0 :                zeta = gto_basis_set%zet(ipgf, iset)
    1040           0 :                gto_basis_set%gcc(ipgf, ishell, iset) = prefac*zeta**expzet*gcca
    1041             :             END DO
    1042           0 :             DO icgf = gto_basis_set%first_cgf(ishell, iset), &
    1043           0 :                gto_basis_set%last_cgf(ishell, iset)
    1044           0 :                gto_basis_set%norm_cgf(icgf) = 1.0_dp
    1045             :             END DO
    1046             :          END DO
    1047             :       END DO
    1048             : 
    1049           0 :    END SUBROUTINE init_norm_cgf_orb_den
    1050             : 
    1051             : ! **************************************************************************************************
    1052             : !> \brief Initialise a Gaussian-type orbital (GTO) basis set data set.
    1053             : !> \param gto_basis_set ...
    1054             : !> \author MK
    1055             : ! **************************************************************************************************
    1056       29822 :    SUBROUTINE init_orb_basis_set(gto_basis_set)
    1057             : 
    1058             :       TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set
    1059             : 
    1060             :       CHARACTER(len=*), PARAMETER :: routineN = 'init_orb_basis_set'
    1061             : 
    1062             :       INTEGER                                            :: handle
    1063             : 
    1064             : ! -------------------------------------------------------------------------
    1065             : 
    1066       14911 :       IF (.NOT. ASSOCIATED(gto_basis_set)) RETURN
    1067             : 
    1068       14911 :       CALL timeset(routineN, handle)
    1069             : 
    1070       14911 :       SELECT CASE (gto_basis_set%norm_type)
    1071             :       CASE (0)
    1072             :          ! No normalisation requested
    1073             :       CASE (1)
    1074           0 :          CALL init_norm_cgf_orb_den(gto_basis_set)
    1075             :       CASE (2)
    1076             :          ! Normalise the primitive Gaussian functions
    1077       14911 :          CALL normalise_gcc_orb(gto_basis_set)
    1078             :          ! Compute the normalization factors of the contracted Gaussian-type
    1079             :          ! functions
    1080       14911 :          CALL init_norm_cgf_orb(gto_basis_set)
    1081             :       CASE DEFAULT
    1082       14911 :          CPABORT("Normalization method not specified")
    1083             :       END SELECT
    1084             : 
    1085             :       ! Initialise the transformation matrices "pgf" -> "cgf"
    1086             : 
    1087       14911 :       CALL init_cphi_and_sphi(gto_basis_set)
    1088             : 
    1089       14911 :       CALL timestop(handle)
    1090             : 
    1091             :    END SUBROUTINE init_orb_basis_set
    1092             : 
    1093             : ! **************************************************************************************************
    1094             : !> \brief Normalise the primitive Cartesian Gaussian functions. The normalization
    1095             : !>        factor is included in the Gaussian contraction coefficients.
    1096             : !> \param gto_basis_set ...
    1097             : !> \author MK
    1098             : ! **************************************************************************************************
    1099       14911 :    SUBROUTINE normalise_gcc_orb(gto_basis_set)
    1100             : 
    1101             :       TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set
    1102             : 
    1103             :       INTEGER                                            :: ipgf, iset, ishell, l
    1104             :       REAL(KIND=dp)                                      :: expzet, gcca, prefac, zeta
    1105             : 
    1106       54496 :       DO iset = 1, gto_basis_set%nset
    1107      120325 :          DO ishell = 1, gto_basis_set%nshell(iset)
    1108       65829 :             l = gto_basis_set%l(ishell, iset)
    1109       65829 :             expzet = 0.25_dp*REAL(2*l + 3, dp)
    1110       65829 :             prefac = 2.0_dp**l*(2.0_dp/pi)**0.75_dp
    1111      352919 :             DO ipgf = 1, gto_basis_set%npgf(iset)
    1112      247505 :                gcca = gto_basis_set%gcc(ipgf, ishell, iset)
    1113      247505 :                zeta = gto_basis_set%zet(ipgf, iset)
    1114      313334 :                gto_basis_set%gcc(ipgf, ishell, iset) = prefac*zeta**expzet*gcca
    1115             :             END DO
    1116             :          END DO
    1117             :       END DO
    1118             : 
    1119       14911 :    END SUBROUTINE normalise_gcc_orb
    1120             : 
    1121             : ! **************************************************************************************************
    1122             : !> \brief Read a Gaussian-type orbital (GTO) basis set from the database file.
    1123             : !> \param element_symbol ...
    1124             : !> \param basis_set_name ...
    1125             : !> \param gto_basis_set ...
    1126             : !> \param para_env ...
    1127             : !> \param dft_section ...
    1128             : !> \author MK
    1129             : ! **************************************************************************************************
    1130       10976 :    SUBROUTINE read_gto_basis_set1(element_symbol, basis_set_name, gto_basis_set, &
    1131             :                                   para_env, dft_section)
    1132             : 
    1133             :       CHARACTER(LEN=*), INTENT(IN)                       :: element_symbol, basis_set_name
    1134             :       TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set
    1135             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1136             :       TYPE(section_vals_type), POINTER                   :: dft_section
    1137             : 
    1138             :       CHARACTER(LEN=240)                                 :: line
    1139             :       CHARACTER(LEN=242)                                 :: line2
    1140             :       CHARACTER(len=default_path_length)                 :: basis_set_file_name, tmp
    1141             :       CHARACTER(LEN=default_path_length), DIMENSION(:), &
    1142       10976 :          POINTER                                         :: cbasis
    1143       10976 :       CHARACTER(LEN=LEN(basis_set_name))                 :: bsname
    1144       10976 :       CHARACTER(LEN=LEN(basis_set_name)+2)               :: bsname2
    1145       10976 :       CHARACTER(LEN=LEN(element_symbol))                 :: symbol
    1146       10976 :       CHARACTER(LEN=LEN(element_symbol)+2)               :: symbol2
    1147             :       INTEGER :: i, ibasis, ico, ipgf, irep, iset, ishell, lshell, m, maxco, maxl, maxpgf, &
    1148             :          maxshell, nbasis, ncgf, nmin, nset, nsgf, sort_method, strlen1, strlen2
    1149       10976 :       INTEGER, DIMENSION(:), POINTER                     :: lmax, lmin, npgf, nshell
    1150       10976 :       INTEGER, DIMENSION(:, :), POINTER                  :: l, n
    1151             :       LOGICAL                                            :: basis_found, found, match
    1152       10976 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zet
    1153       10976 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: gcc
    1154             :       TYPE(cp_parser_type)                               :: parser
    1155             : 
    1156       10976 :       line = ""
    1157       10976 :       line2 = ""
    1158       10976 :       symbol = ""
    1159       10976 :       symbol2 = ""
    1160       10976 :       bsname = ""
    1161       10976 :       bsname2 = ""
    1162             : 
    1163       10976 :       nbasis = 1
    1164             : 
    1165       10976 :       gto_basis_set%name = basis_set_name
    1166       10976 :       gto_basis_set%aliases = basis_set_name
    1167             :       CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", &
    1168       10976 :                                 n_rep_val=nbasis)
    1169       32928 :       ALLOCATE (cbasis(nbasis))
    1170       23804 :       DO ibasis = 1, nbasis
    1171             :          CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", &
    1172       12828 :                                    i_rep_val=ibasis, c_val=cbasis(ibasis))
    1173       12828 :          basis_set_file_name = cbasis(ibasis)
    1174       12828 :          tmp = basis_set_file_name
    1175       12828 :          CALL uppercase(tmp)
    1176       23804 :          IF (INDEX(tmp, "MOLOPT") .NE. 0) CALL cite_reference(VandeVondele2007)
    1177             :       END DO
    1178             : 
    1179             :       ! Search for the requested basis set in the basis set file
    1180             :       ! until the basis set is found or the end of file is reached
    1181             : 
    1182       10976 :       basis_found = .FALSE.
    1183       22804 :       basis_loop: DO ibasis = 1, nbasis
    1184       12758 :          IF (basis_found) EXIT basis_loop
    1185       11828 :          basis_set_file_name = cbasis(ibasis)
    1186       11828 :          CALL parser_create(parser, basis_set_file_name, para_env=para_env)
    1187             : 
    1188       11828 :          bsname = basis_set_name
    1189       11828 :          symbol = element_symbol
    1190       11828 :          irep = 0
    1191             : 
    1192       11828 :          tmp = basis_set_name
    1193       11828 :          CALL uppercase(tmp)
    1194       11828 :          IF (INDEX(tmp, "MOLOPT") .NE. 0) CALL cite_reference(VandeVondele2007)
    1195             : 
    1196       11828 :          nset = 0
    1197       11828 :          maxshell = 0
    1198       11828 :          maxpgf = 0
    1199       11828 :          maxco = 0
    1200       11828 :          ncgf = 0
    1201       11828 :          nsgf = 0
    1202       11828 :          gto_basis_set%nset = nset
    1203       11828 :          gto_basis_set%ncgf = ncgf
    1204       11828 :          gto_basis_set%nsgf = nsgf
    1205       11828 :          CALL reallocate(gto_basis_set%lmax, 1, nset)
    1206       11828 :          CALL reallocate(gto_basis_set%lmin, 1, nset)
    1207       11828 :          CALL reallocate(gto_basis_set%npgf, 1, nset)
    1208       11828 :          CALL reallocate(gto_basis_set%nshell, 1, nset)
    1209       11828 :          CALL reallocate(gto_basis_set%n, 1, maxshell, 1, nset)
    1210       11828 :          CALL reallocate(gto_basis_set%l, 1, maxshell, 1, nset)
    1211       11828 :          CALL reallocate(gto_basis_set%zet, 1, maxpgf, 1, nset)
    1212       11828 :          CALL reallocate(gto_basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)
    1213       11828 :          CALL reallocate(gto_basis_set%set_radius, 1, nset)
    1214       11828 :          CALL reallocate(gto_basis_set%pgf_radius, 1, maxpgf, 1, nset)
    1215       11828 :          CALL reallocate(gto_basis_set%first_cgf, 1, maxshell, 1, nset)
    1216       11828 :          CALL reallocate(gto_basis_set%first_sgf, 1, maxshell, 1, nset)
    1217       11828 :          CALL reallocate(gto_basis_set%last_cgf, 1, maxshell, 1, nset)
    1218       11828 :          CALL reallocate(gto_basis_set%last_sgf, 1, maxshell, 1, nset)
    1219       11828 :          CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
    1220       11828 :          CALL reallocate(gto_basis_set%nsgf_set, 1, nset)
    1221       11828 :          CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
    1222       11828 :          CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
    1223       11828 :          CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
    1224       11828 :          CALL reallocate(gto_basis_set%lx, 1, ncgf)
    1225       11828 :          CALL reallocate(gto_basis_set%ly, 1, ncgf)
    1226       11828 :          CALL reallocate(gto_basis_set%lz, 1, ncgf)
    1227       11828 :          CALL reallocate(gto_basis_set%m, 1, nsgf)
    1228       11828 :          CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
    1229             : 
    1230       11828 :          IF (tmp .NE. "NONE") THEN
    1231             :             search_loop: DO
    1232             : 
    1233       59057 :                CALL parser_search_string(parser, TRIM(bsname), .TRUE., found, line)
    1234       59057 :                IF (found) THEN
    1235       58205 :                   CALL uppercase(symbol)
    1236       58205 :                   CALL uppercase(bsname)
    1237             : 
    1238       58205 :                   match = .FALSE.
    1239       58205 :                   CALL uppercase(line)
    1240             :                   ! Check both the element symbol and the basis set name
    1241       58205 :                   line2 = " "//line//" "
    1242       58205 :                   symbol2 = " "//TRIM(symbol)//" "
    1243       58205 :                   bsname2 = " "//TRIM(bsname)//" "
    1244       58205 :                   strlen1 = LEN_TRIM(symbol2) + 1
    1245       58205 :                   strlen2 = LEN_TRIM(bsname2) + 1
    1246             : 
    1247       58205 :                   IF ((INDEX(line2, symbol2(:strlen1)) > 0) .AND. &
    1248             :                       (INDEX(line2, bsname2(:strlen2)) > 0)) match = .TRUE.
    1249             :                   IF (match) THEN
    1250             :                      ! copy all names into aliases field
    1251       10968 :                      i = INDEX(line2, symbol2(:strlen1))
    1252       10968 :                      i = i + 1 + INDEX(line2(i + 1:), " ")
    1253       10968 :                      gto_basis_set%aliases = line2(i:)
    1254             : 
    1255       10968 :                      NULLIFY (gcc, l, lmax, lmin, n, npgf, nshell, zet)
    1256             :                      ! Read the basis set information
    1257       10968 :                      CALL parser_get_object(parser, nset, newline=.TRUE.)
    1258             : 
    1259       10968 :                      CALL reallocate(npgf, 1, nset)
    1260       10968 :                      CALL reallocate(nshell, 1, nset)
    1261       10968 :                      CALL reallocate(lmax, 1, nset)
    1262       10968 :                      CALL reallocate(lmin, 1, nset)
    1263       10968 :                      CALL reallocate(n, 1, 1, 1, nset)
    1264             : 
    1265       10968 :                      maxl = 0
    1266       10968 :                      maxpgf = 0
    1267       10968 :                      maxshell = 0
    1268             : 
    1269       42360 :                      DO iset = 1, nset
    1270       31392 :                         CALL parser_get_object(parser, n(1, iset), newline=.TRUE.)
    1271       31392 :                         CALL parser_get_object(parser, lmin(iset))
    1272       31392 :                         CALL parser_get_object(parser, lmax(iset))
    1273       31392 :                         CALL parser_get_object(parser, npgf(iset))
    1274       31392 :                         maxl = MAX(maxl, lmax(iset))
    1275       31392 :                         IF (npgf(iset) > maxpgf) THEN
    1276       11156 :                            maxpgf = npgf(iset)
    1277       11156 :                            CALL reallocate(zet, 1, maxpgf, 1, nset)
    1278       11156 :                            CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
    1279             :                         END IF
    1280       31392 :                         nshell(iset) = 0
    1281       71526 :                         DO lshell = lmin(iset), lmax(iset)
    1282       40134 :                            nmin = n(1, iset) + lshell - lmin(iset)
    1283       40134 :                            CALL parser_get_object(parser, ishell)
    1284       40134 :                            nshell(iset) = nshell(iset) + ishell
    1285       40134 :                            IF (nshell(iset) > maxshell) THEN
    1286       17142 :                               maxshell = nshell(iset)
    1287       17142 :                               CALL reallocate(n, 1, maxshell, 1, nset)
    1288       17142 :                               CALL reallocate(l, 1, maxshell, 1, nset)
    1289       17142 :                               CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
    1290             :                            END IF
    1291      162193 :                            DO i = 1, ishell
    1292       50533 :                               n(nshell(iset) - ishell + i, iset) = nmin + i - 1
    1293       90667 :                               l(nshell(iset) - ishell + i, iset) = lshell
    1294             :                            END DO
    1295             :                         END DO
    1296      108972 :                         DO ipgf = 1, npgf(iset)
    1297       66612 :                            CALL parser_get_object(parser, zet(ipgf, iset), newline=.TRUE.)
    1298      243579 :                            DO ishell = 1, nshell(iset)
    1299      212187 :                               CALL parser_get_object(parser, gcc(ipgf, ishell, iset))
    1300             :                            END DO
    1301             :                         END DO
    1302             :                      END DO
    1303             : 
    1304             :                      ! Maximum angular momentum quantum number of the atomic kind
    1305             : 
    1306       10968 :                      CALL init_orbital_pointers(maxl)
    1307             : 
    1308             :                      ! Allocate the global variables
    1309             : 
    1310       10968 :                      gto_basis_set%nset = nset
    1311       10968 :                      CALL reallocate(gto_basis_set%lmax, 1, nset)
    1312       10968 :                      CALL reallocate(gto_basis_set%lmin, 1, nset)
    1313       10968 :                      CALL reallocate(gto_basis_set%npgf, 1, nset)
    1314       10968 :                      CALL reallocate(gto_basis_set%nshell, 1, nset)
    1315       10968 :                      CALL reallocate(gto_basis_set%n, 1, maxshell, 1, nset)
    1316       10968 :                      CALL reallocate(gto_basis_set%l, 1, maxshell, 1, nset)
    1317       10968 :                      CALL reallocate(gto_basis_set%zet, 1, maxpgf, 1, nset)
    1318       10968 :                      CALL reallocate(gto_basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)
    1319             : 
    1320             :                      ! Copy the basis set information into the data structure
    1321             : 
    1322       42360 :                      DO iset = 1, nset
    1323       31392 :                         gto_basis_set%lmax(iset) = lmax(iset)
    1324       31392 :                         gto_basis_set%lmin(iset) = lmin(iset)
    1325       31392 :                         gto_basis_set%npgf(iset) = npgf(iset)
    1326       31392 :                         gto_basis_set%nshell(iset) = nshell(iset)
    1327       81925 :                         DO ishell = 1, nshell(iset)
    1328       50533 :                            gto_basis_set%n(ishell, iset) = n(ishell, iset)
    1329       50533 :                            gto_basis_set%l(ishell, iset) = l(ishell, iset)
    1330      227500 :                            DO ipgf = 1, npgf(iset)
    1331      196108 :                               gto_basis_set%gcc(ipgf, ishell, iset) = gcc(ipgf, ishell, iset)
    1332             :                            END DO
    1333             :                         END DO
    1334      108972 :                         DO ipgf = 1, npgf(iset)
    1335       98004 :                            gto_basis_set%zet(ipgf, iset) = zet(ipgf, iset)
    1336             :                         END DO
    1337             :                      END DO
    1338             : 
    1339             :                      ! Initialise the depending atomic kind information
    1340             : 
    1341       10968 :                      CALL reallocate(gto_basis_set%set_radius, 1, nset)
    1342       10968 :                      CALL reallocate(gto_basis_set%pgf_radius, 1, maxpgf, 1, nset)
    1343       10968 :                      CALL reallocate(gto_basis_set%first_cgf, 1, maxshell, 1, nset)
    1344       10968 :                      CALL reallocate(gto_basis_set%first_sgf, 1, maxshell, 1, nset)
    1345       10968 :                      CALL reallocate(gto_basis_set%last_cgf, 1, maxshell, 1, nset)
    1346       10968 :                      CALL reallocate(gto_basis_set%last_sgf, 1, maxshell, 1, nset)
    1347       10968 :                      CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
    1348       10968 :                      CALL reallocate(gto_basis_set%nsgf_set, 1, nset)
    1349             : 
    1350       10968 :                      maxco = 0
    1351       10968 :                      ncgf = 0
    1352       10968 :                      nsgf = 0
    1353             : 
    1354       42360 :                      DO iset = 1, nset
    1355       31392 :                         gto_basis_set%ncgf_set(iset) = 0
    1356       31392 :                         gto_basis_set%nsgf_set(iset) = 0
    1357       81925 :                         DO ishell = 1, nshell(iset)
    1358       50533 :                            lshell = gto_basis_set%l(ishell, iset)
    1359       50533 :                            gto_basis_set%first_cgf(ishell, iset) = ncgf + 1
    1360       50533 :                            ncgf = ncgf + nco(lshell)
    1361       50533 :                            gto_basis_set%last_cgf(ishell, iset) = ncgf
    1362             :                            gto_basis_set%ncgf_set(iset) = &
    1363       50533 :                               gto_basis_set%ncgf_set(iset) + nco(lshell)
    1364       50533 :                            gto_basis_set%first_sgf(ishell, iset) = nsgf + 1
    1365       50533 :                            nsgf = nsgf + nso(lshell)
    1366       50533 :                            gto_basis_set%last_sgf(ishell, iset) = nsgf
    1367             :                            gto_basis_set%nsgf_set(iset) = &
    1368       81925 :                               gto_basis_set%nsgf_set(iset) + nso(lshell)
    1369             :                         END DO
    1370       42360 :                         maxco = MAX(maxco, npgf(iset)*ncoset(lmax(iset)))
    1371             :                      END DO
    1372             : 
    1373       10968 :                      gto_basis_set%ncgf = ncgf
    1374       10968 :                      gto_basis_set%nsgf = nsgf
    1375             : 
    1376       10968 :                      CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
    1377       10968 :                      CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
    1378       10968 :                      CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
    1379       10968 :                      CALL reallocate(gto_basis_set%lx, 1, ncgf)
    1380       10968 :                      CALL reallocate(gto_basis_set%ly, 1, ncgf)
    1381       10968 :                      CALL reallocate(gto_basis_set%lz, 1, ncgf)
    1382       10968 :                      CALL reallocate(gto_basis_set%m, 1, nsgf)
    1383       10968 :                      CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
    1384       32904 :                      ALLOCATE (gto_basis_set%cgf_symbol(ncgf))
    1385             : 
    1386       32904 :                      ALLOCATE (gto_basis_set%sgf_symbol(nsgf))
    1387             : 
    1388       10968 :                      ncgf = 0
    1389       10968 :                      nsgf = 0
    1390             : 
    1391       42360 :                      DO iset = 1, nset
    1392       92893 :                         DO ishell = 1, nshell(iset)
    1393       50533 :                            lshell = gto_basis_set%l(ishell, iset)
    1394      188183 :                            DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
    1395      137650 :                               ncgf = ncgf + 1
    1396      137650 :                               gto_basis_set%lx(ncgf) = indco(1, ico)
    1397      137650 :                               gto_basis_set%ly(ncgf) = indco(2, ico)
    1398      137650 :                               gto_basis_set%lz(ncgf) = indco(3, ico)
    1399             :                               gto_basis_set%cgf_symbol(ncgf) = &
    1400             :                                  cgf_symbol(n(ishell, iset), (/gto_basis_set%lx(ncgf), &
    1401             :                                                                gto_basis_set%ly(ncgf), &
    1402      601133 :                                                                gto_basis_set%lz(ncgf)/))
    1403             :                            END DO
    1404      206320 :                            DO m = -lshell, lshell
    1405      124395 :                               nsgf = nsgf + 1
    1406      124395 :                               gto_basis_set%m(nsgf) = m
    1407             :                               gto_basis_set%sgf_symbol(nsgf) = &
    1408      174928 :                                  sgf_symbol(n(ishell, iset), lshell, m)
    1409             :                            END DO
    1410             :                         END DO
    1411             :                      END DO
    1412             : 
    1413       10968 :                      DEALLOCATE (gcc, l, lmax, lmin, n, npgf, nshell, zet)
    1414             : 
    1415       10968 :                      basis_found = .TRUE.
    1416       10968 :                      EXIT search_loop
    1417             :                   END IF
    1418             :                ELSE
    1419             :                   EXIT search_loop
    1420             :                END IF
    1421             :             END DO search_loop
    1422             :          ELSE
    1423           8 :             match = .FALSE.
    1424          16 :             ALLOCATE (gto_basis_set%cgf_symbol(ncgf))
    1425          16 :             ALLOCATE (gto_basis_set%sgf_symbol(nsgf))
    1426             :          END IF
    1427             : 
    1428       22804 :          CALL parser_release(parser)
    1429             : 
    1430             :       END DO basis_loop
    1431             : 
    1432       10976 :       IF (tmp .NE. "NONE") THEN
    1433       10968 :          IF (.NOT. basis_found) THEN
    1434           0 :             basis_set_file_name = ""
    1435           0 :             DO ibasis = 1, nbasis
    1436           0 :                basis_set_file_name = TRIM(basis_set_file_name)//"<"//TRIM(cbasis(ibasis))//"> "
    1437             :             END DO
    1438             :             CALL cp_abort(__LOCATION__, &
    1439             :                           "The requested basis set <"//TRIM(bsname)// &
    1440             :                           "> for element <"//TRIM(symbol)//"> was not "// &
    1441             :                           "found in the basis set files "// &
    1442           0 :                           TRIM(basis_set_file_name))
    1443             :          END IF
    1444             :       END IF
    1445       10976 :       DEALLOCATE (cbasis)
    1446             : 
    1447       10976 :       CALL section_vals_val_get(dft_section, "SORT_BASIS", i_val=sort_method)
    1448       10976 :       CALL sort_gto_basis_set(gto_basis_set, sort_method)
    1449             : 
    1450       54880 :    END SUBROUTINE read_gto_basis_set1
    1451             : 
    1452             : ! **************************************************************************************************
    1453             : !> \brief Read a Gaussian-type orbital (GTO) basis set from the database file.
    1454             : !> \param element_symbol ...
    1455             : !> \param basis_type ...
    1456             : !> \param gto_basis_set ...
    1457             : !> \param basis_section ...
    1458             : !> \param irep ...
    1459             : !> \param dft_section ...
    1460             : !> \author MK
    1461             : ! **************************************************************************************************
    1462         174 :    SUBROUTINE read_gto_basis_set2(element_symbol, basis_type, gto_basis_set, &
    1463             :                                   basis_section, irep, dft_section)
    1464             : 
    1465             :       CHARACTER(LEN=*), INTENT(IN)                       :: element_symbol
    1466             :       CHARACTER(LEN=*), INTENT(INOUT)                    :: basis_type
    1467             :       TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set
    1468             :       TYPE(section_vals_type), OPTIONAL, POINTER         :: basis_section
    1469             :       INTEGER                                            :: irep
    1470             :       TYPE(section_vals_type), OPTIONAL, POINTER         :: dft_section
    1471             : 
    1472             :       CHARACTER(len=20*default_string_length)            :: line_att
    1473             :       CHARACTER(LEN=240)                                 :: line
    1474             :       CHARACTER(LEN=242)                                 :: line2
    1475             :       CHARACTER(LEN=default_path_length)                 :: bsname, bsname2
    1476         174 :       CHARACTER(LEN=LEN(element_symbol))                 :: symbol
    1477         174 :       CHARACTER(LEN=LEN(element_symbol)+2)               :: symbol2
    1478             :       INTEGER                                            :: i, ico, ipgf, iset, ishell, lshell, m, &
    1479             :                                                             maxco, maxl, maxpgf, maxshell, ncgf, &
    1480             :                                                             nmin, nset, nsgf, sort_method
    1481         174 :       INTEGER, DIMENSION(:), POINTER                     :: lmax, lmin, npgf, nshell
    1482         174 :       INTEGER, DIMENSION(:, :), POINTER                  :: l, n
    1483             :       LOGICAL                                            :: is_ok
    1484         174 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: zet
    1485         174 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: gcc
    1486             :       TYPE(cp_sll_val_type), POINTER                     :: list
    1487             :       TYPE(val_type), POINTER                            :: val
    1488             : 
    1489         174 :       line = ""
    1490         174 :       line2 = ""
    1491         174 :       symbol = ""
    1492         174 :       symbol2 = ""
    1493             :       bsname = ""
    1494         174 :       bsname2 = ""
    1495             : 
    1496         174 :       gto_basis_set%name = " "
    1497         174 :       gto_basis_set%aliases = " "
    1498             : 
    1499         174 :       bsname = " "
    1500         174 :       symbol = element_symbol
    1501             : 
    1502         174 :       nset = 0
    1503         174 :       maxshell = 0
    1504         174 :       maxpgf = 0
    1505         174 :       maxco = 0
    1506         174 :       ncgf = 0
    1507         174 :       nsgf = 0
    1508         174 :       gto_basis_set%nset = nset
    1509         174 :       gto_basis_set%ncgf = ncgf
    1510         174 :       gto_basis_set%nsgf = nsgf
    1511         174 :       CALL reallocate(gto_basis_set%lmax, 1, nset)
    1512         174 :       CALL reallocate(gto_basis_set%lmin, 1, nset)
    1513         174 :       CALL reallocate(gto_basis_set%npgf, 1, nset)
    1514         174 :       CALL reallocate(gto_basis_set%nshell, 1, nset)
    1515         174 :       CALL reallocate(gto_basis_set%n, 1, maxshell, 1, nset)
    1516         174 :       CALL reallocate(gto_basis_set%l, 1, maxshell, 1, nset)
    1517         174 :       CALL reallocate(gto_basis_set%zet, 1, maxpgf, 1, nset)
    1518         174 :       CALL reallocate(gto_basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)
    1519         174 :       CALL reallocate(gto_basis_set%set_radius, 1, nset)
    1520         174 :       CALL reallocate(gto_basis_set%pgf_radius, 1, maxpgf, 1, nset)
    1521         174 :       CALL reallocate(gto_basis_set%first_cgf, 1, maxshell, 1, nset)
    1522         174 :       CALL reallocate(gto_basis_set%first_sgf, 1, maxshell, 1, nset)
    1523         174 :       CALL reallocate(gto_basis_set%last_cgf, 1, maxshell, 1, nset)
    1524         174 :       CALL reallocate(gto_basis_set%last_sgf, 1, maxshell, 1, nset)
    1525         174 :       CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
    1526         174 :       CALL reallocate(gto_basis_set%nsgf_set, 1, nset)
    1527         174 :       CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
    1528         174 :       CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
    1529         174 :       CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
    1530         174 :       CALL reallocate(gto_basis_set%lx, 1, ncgf)
    1531         174 :       CALL reallocate(gto_basis_set%ly, 1, ncgf)
    1532         174 :       CALL reallocate(gto_basis_set%lz, 1, ncgf)
    1533         174 :       CALL reallocate(gto_basis_set%m, 1, nsgf)
    1534         174 :       CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
    1535             : 
    1536         174 :       basis_type = ""
    1537         174 :       CALL section_vals_val_get(basis_section, "_SECTION_PARAMETERS_", i_rep_section=irep, c_val=basis_type)
    1538         174 :       IF (basis_type == "Orbital") basis_type = "ORB"
    1539             : 
    1540         174 :       NULLIFY (list, val)
    1541         174 :       CALL section_vals_list_get(basis_section, "_DEFAULT_KEYWORD_", i_rep_section=irep, list=list)
    1542         174 :       CALL uppercase(symbol)
    1543         174 :       CALL uppercase(bsname)
    1544             : 
    1545         174 :       NULLIFY (gcc, l, lmax, lmin, n, npgf, nshell, zet)
    1546             :       ! Read the basis set information
    1547         174 :       is_ok = cp_sll_val_next(list, val)
    1548         174 :       IF (.NOT. is_ok) CPABORT("Error reading the Basis set from input file!!")
    1549         174 :       CALL val_get(val, c_val=line_att)
    1550         174 :       READ (line_att, *) nset
    1551             : 
    1552         174 :       CALL reallocate(npgf, 1, nset)
    1553         174 :       CALL reallocate(nshell, 1, nset)
    1554         174 :       CALL reallocate(lmax, 1, nset)
    1555         174 :       CALL reallocate(lmin, 1, nset)
    1556         174 :       CALL reallocate(n, 1, 1, 1, nset)
    1557             : 
    1558         174 :       maxl = 0
    1559         174 :       maxpgf = 0
    1560         174 :       maxshell = 0
    1561             : 
    1562         500 :       DO iset = 1, nset
    1563         326 :          is_ok = cp_sll_val_next(list, val)
    1564         326 :          IF (.NOT. is_ok) CPABORT("Error reading the Basis set from input file!!")
    1565         326 :          CALL val_get(val, c_val=line_att)
    1566         326 :          READ (line_att, *) n(1, iset)
    1567         326 :          CALL remove_word(line_att)
    1568         326 :          READ (line_att, *) lmin(iset)
    1569         326 :          CALL remove_word(line_att)
    1570         326 :          READ (line_att, *) lmax(iset)
    1571         326 :          CALL remove_word(line_att)
    1572         326 :          READ (line_att, *) npgf(iset)
    1573         326 :          CALL remove_word(line_att)
    1574         326 :          maxl = MAX(maxl, lmax(iset))
    1575         326 :          IF (npgf(iset) > maxpgf) THEN
    1576         174 :             maxpgf = npgf(iset)
    1577         174 :             CALL reallocate(zet, 1, maxpgf, 1, nset)
    1578         174 :             CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
    1579             :          END IF
    1580         326 :          nshell(iset) = 0
    1581        1082 :          DO lshell = lmin(iset), lmax(iset)
    1582         756 :             nmin = n(1, iset) + lshell - lmin(iset)
    1583         756 :             READ (line_att, *) ishell
    1584         756 :             CALL remove_word(line_att)
    1585         756 :             nshell(iset) = nshell(iset) + ishell
    1586         756 :             IF (nshell(iset) > maxshell) THEN
    1587         334 :                maxshell = nshell(iset)
    1588         334 :                CALL reallocate(n, 1, maxshell, 1, nset)
    1589         334 :                CALL reallocate(l, 1, maxshell, 1, nset)
    1590         334 :                CALL reallocate(gcc, 1, maxpgf, 1, maxshell, 1, nset)
    1591             :             END IF
    1592        2126 :             DO i = 1, ishell
    1593        1044 :                n(nshell(iset) - ishell + i, iset) = nmin + i - 1
    1594        1800 :                l(nshell(iset) - ishell + i, iset) = lshell
    1595             :             END DO
    1596             :          END DO
    1597         326 :          IF (LEN_TRIM(line_att) /= 0) &
    1598           0 :             CPABORT("Error reading the Basis from input file!!")
    1599        1308 :          DO ipgf = 1, npgf(iset)
    1600         808 :             is_ok = cp_sll_val_next(list, val)
    1601         808 :             IF (.NOT. is_ok) CPABORT("Error reading the Basis set from input file!!")
    1602         808 :             CALL val_get(val, c_val=line_att)
    1603        1134 :             READ (line_att, *) zet(ipgf, iset), (gcc(ipgf, ishell, iset), ishell=1, nshell(iset))
    1604             :          END DO
    1605             :       END DO
    1606             : 
    1607             :       ! Maximum angular momentum quantum number of the atomic kind
    1608             : 
    1609         174 :       CALL init_orbital_pointers(maxl)
    1610             : 
    1611             :       ! Allocate the global variables
    1612             : 
    1613         174 :       gto_basis_set%nset = nset
    1614         174 :       CALL reallocate(gto_basis_set%lmax, 1, nset)
    1615         174 :       CALL reallocate(gto_basis_set%lmin, 1, nset)
    1616         174 :       CALL reallocate(gto_basis_set%npgf, 1, nset)
    1617         174 :       CALL reallocate(gto_basis_set%nshell, 1, nset)
    1618         174 :       CALL reallocate(gto_basis_set%n, 1, maxshell, 1, nset)
    1619         174 :       CALL reallocate(gto_basis_set%l, 1, maxshell, 1, nset)
    1620         174 :       CALL reallocate(gto_basis_set%zet, 1, maxpgf, 1, nset)
    1621         174 :       CALL reallocate(gto_basis_set%gcc, 1, maxpgf, 1, maxshell, 1, nset)
    1622             : 
    1623             :       ! Copy the basis set information into the data structure
    1624             : 
    1625         500 :       DO iset = 1, nset
    1626         326 :          gto_basis_set%lmax(iset) = lmax(iset)
    1627         326 :          gto_basis_set%lmin(iset) = lmin(iset)
    1628         326 :          gto_basis_set%npgf(iset) = npgf(iset)
    1629         326 :          gto_basis_set%nshell(iset) = nshell(iset)
    1630        1370 :          DO ishell = 1, nshell(iset)
    1631        1044 :             gto_basis_set%n(ishell, iset) = n(ishell, iset)
    1632        1044 :             gto_basis_set%l(ishell, iset) = l(ishell, iset)
    1633        5656 :             DO ipgf = 1, npgf(iset)
    1634        5330 :                gto_basis_set%gcc(ipgf, ishell, iset) = gcc(ipgf, ishell, iset)
    1635             :             END DO
    1636             :          END DO
    1637        1308 :          DO ipgf = 1, npgf(iset)
    1638        1134 :             gto_basis_set%zet(ipgf, iset) = zet(ipgf, iset)
    1639             :          END DO
    1640             :       END DO
    1641             : 
    1642             :       ! Initialise the depending atomic kind information
    1643             : 
    1644         174 :       CALL reallocate(gto_basis_set%set_radius, 1, nset)
    1645         174 :       CALL reallocate(gto_basis_set%pgf_radius, 1, maxpgf, 1, nset)
    1646         174 :       CALL reallocate(gto_basis_set%first_cgf, 1, maxshell, 1, nset)
    1647         174 :       CALL reallocate(gto_basis_set%first_sgf, 1, maxshell, 1, nset)
    1648         174 :       CALL reallocate(gto_basis_set%last_cgf, 1, maxshell, 1, nset)
    1649         174 :       CALL reallocate(gto_basis_set%last_sgf, 1, maxshell, 1, nset)
    1650         174 :       CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
    1651         174 :       CALL reallocate(gto_basis_set%nsgf_set, 1, nset)
    1652             : 
    1653         174 :       maxco = 0
    1654         174 :       ncgf = 0
    1655         174 :       nsgf = 0
    1656             : 
    1657         500 :       DO iset = 1, nset
    1658         326 :          gto_basis_set%ncgf_set(iset) = 0
    1659         326 :          gto_basis_set%nsgf_set(iset) = 0
    1660        1370 :          DO ishell = 1, nshell(iset)
    1661        1044 :             lshell = gto_basis_set%l(ishell, iset)
    1662        1044 :             gto_basis_set%first_cgf(ishell, iset) = ncgf + 1
    1663        1044 :             ncgf = ncgf + nco(lshell)
    1664        1044 :             gto_basis_set%last_cgf(ishell, iset) = ncgf
    1665             :             gto_basis_set%ncgf_set(iset) = &
    1666        1044 :                gto_basis_set%ncgf_set(iset) + nco(lshell)
    1667        1044 :             gto_basis_set%first_sgf(ishell, iset) = nsgf + 1
    1668        1044 :             nsgf = nsgf + nso(lshell)
    1669        1044 :             gto_basis_set%last_sgf(ishell, iset) = nsgf
    1670             :             gto_basis_set%nsgf_set(iset) = &
    1671        1370 :                gto_basis_set%nsgf_set(iset) + nso(lshell)
    1672             :          END DO
    1673         500 :          maxco = MAX(maxco, npgf(iset)*ncoset(lmax(iset)))
    1674             :       END DO
    1675             : 
    1676         174 :       gto_basis_set%ncgf = ncgf
    1677         174 :       gto_basis_set%nsgf = nsgf
    1678             : 
    1679         174 :       CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
    1680         174 :       CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
    1681         174 :       CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
    1682         174 :       CALL reallocate(gto_basis_set%lx, 1, ncgf)
    1683         174 :       CALL reallocate(gto_basis_set%ly, 1, ncgf)
    1684         174 :       CALL reallocate(gto_basis_set%lz, 1, ncgf)
    1685         174 :       CALL reallocate(gto_basis_set%m, 1, nsgf)
    1686         174 :       CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
    1687         522 :       ALLOCATE (gto_basis_set%cgf_symbol(ncgf))
    1688             : 
    1689         522 :       ALLOCATE (gto_basis_set%sgf_symbol(nsgf))
    1690             : 
    1691         174 :       ncgf = 0
    1692         174 :       nsgf = 0
    1693             : 
    1694         500 :       DO iset = 1, nset
    1695        1544 :          DO ishell = 1, nshell(iset)
    1696        1044 :             lshell = gto_basis_set%l(ishell, iset)
    1697        5040 :             DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
    1698        3996 :                ncgf = ncgf + 1
    1699        3996 :                gto_basis_set%lx(ncgf) = indco(1, ico)
    1700        3996 :                gto_basis_set%ly(ncgf) = indco(2, ico)
    1701        3996 :                gto_basis_set%lz(ncgf) = indco(3, ico)
    1702             :                gto_basis_set%cgf_symbol(ncgf) = &
    1703             :                   cgf_symbol(n(ishell, iset), (/gto_basis_set%lx(ncgf), &
    1704             :                                                 gto_basis_set%ly(ncgf), &
    1705       17028 :                                                 gto_basis_set%lz(ncgf)/))
    1706             :             END DO
    1707        4606 :             DO m = -lshell, lshell
    1708        3236 :                nsgf = nsgf + 1
    1709        3236 :                gto_basis_set%m(nsgf) = m
    1710             :                gto_basis_set%sgf_symbol(nsgf) = &
    1711        4280 :                   sgf_symbol(n(ishell, iset), lshell, m)
    1712             :             END DO
    1713             :          END DO
    1714             :       END DO
    1715             : 
    1716         174 :       DEALLOCATE (gcc, l, lmax, lmin, n, npgf, nshell, zet)
    1717             : 
    1718         174 :       IF (PRESENT(dft_section)) THEN
    1719         162 :          CALL section_vals_val_get(dft_section, "SORT_BASIS", i_val=sort_method)
    1720         162 :          CALL sort_gto_basis_set(gto_basis_set, sort_method)
    1721             :       END IF
    1722             : 
    1723         174 :    END SUBROUTINE read_gto_basis_set2
    1724             : 
    1725             : ! **************************************************************************************************
    1726             : !> \brief Set the components of Gaussian-type orbital (GTO) basis set data set.
    1727             : !> \param gto_basis_set ...
    1728             : !> \param name ...
    1729             : !> \param aliases ...
    1730             : !> \param norm_type ...
    1731             : !> \param kind_radius ...
    1732             : !> \param ncgf ...
    1733             : !> \param nset ...
    1734             : !> \param nsgf ...
    1735             : !> \param cgf_symbol ...
    1736             : !> \param sgf_symbol ...
    1737             : !> \param norm_cgf ...
    1738             : !> \param set_radius ...
    1739             : !> \param lmax ...
    1740             : !> \param lmin ...
    1741             : !> \param lx ...
    1742             : !> \param ly ...
    1743             : !> \param lz ...
    1744             : !> \param m ...
    1745             : !> \param ncgf_set ...
    1746             : !> \param npgf ...
    1747             : !> \param nsgf_set ...
    1748             : !> \param nshell ...
    1749             : !> \param cphi ...
    1750             : !> \param pgf_radius ...
    1751             : !> \param sphi ...
    1752             : !> \param scon ...
    1753             : !> \param zet ...
    1754             : !> \param first_cgf ...
    1755             : !> \param first_sgf ...
    1756             : !> \param l ...
    1757             : !> \param last_cgf ...
    1758             : !> \param last_sgf ...
    1759             : !> \param n ...
    1760             : !> \param gcc ...
    1761             : !> \param short_kind_radius ...
    1762             : !> \author MK
    1763             : ! **************************************************************************************************
    1764       31992 :    SUBROUTINE set_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, &
    1765             :                                 nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, &
    1766             :                                 lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, &
    1767             :                                 cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, &
    1768             :                                 last_cgf, last_sgf, n, gcc, short_kind_radius)
    1769             : 
    1770             :       TYPE(gto_basis_set_type), INTENT(INOUT)            :: gto_basis_set
    1771             :       CHARACTER(LEN=default_string_length), INTENT(IN), &
    1772             :          OPTIONAL                                        :: name, aliases
    1773             :       INTEGER, INTENT(IN), OPTIONAL                      :: norm_type
    1774             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: kind_radius
    1775             :       INTEGER, INTENT(IN), OPTIONAL                      :: ncgf, nset, nsgf
    1776             :       CHARACTER(LEN=12), DIMENSION(:), OPTIONAL, POINTER :: cgf_symbol
    1777             :       CHARACTER(LEN=6), DIMENSION(:), OPTIONAL, POINTER  :: sgf_symbol
    1778             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: norm_cgf, set_radius
    1779             :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: lmax, lmin, lx, ly, lz, m, ncgf_set, &
    1780             :                                                             npgf, nsgf_set, nshell
    1781             :       REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: cphi, pgf_radius, sphi, scon, zet
    1782             :       INTEGER, DIMENSION(:, :), OPTIONAL, POINTER        :: first_cgf, first_sgf, l, last_cgf, &
    1783             :                                                             last_sgf, n
    1784             :       REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
    1785             :          POINTER                                         :: gcc
    1786             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: short_kind_radius
    1787             : 
    1788       31992 :       IF (PRESENT(name)) gto_basis_set%name = name
    1789       31992 :       IF (PRESENT(aliases)) gto_basis_set%aliases = aliases
    1790       31992 :       IF (PRESENT(norm_type)) gto_basis_set%norm_type = norm_type
    1791       31992 :       IF (PRESENT(kind_radius)) gto_basis_set%kind_radius = kind_radius
    1792       31992 :       IF (PRESENT(short_kind_radius)) gto_basis_set%short_kind_radius = short_kind_radius
    1793       31992 :       IF (PRESENT(ncgf)) gto_basis_set%ncgf = ncgf
    1794       31992 :       IF (PRESENT(nset)) gto_basis_set%nset = nset
    1795       31992 :       IF (PRESENT(nsgf)) gto_basis_set%nsgf = nsgf
    1796       31992 :       IF (PRESENT(cgf_symbol)) gto_basis_set%cgf_symbol(:) = cgf_symbol(:)
    1797       31992 :       IF (PRESENT(sgf_symbol)) gto_basis_set%sgf_symbol(:) = sgf_symbol(:)
    1798       31992 :       IF (PRESENT(norm_cgf)) gto_basis_set%norm_cgf(:) = norm_cgf(:)
    1799      143270 :       IF (PRESENT(set_radius)) gto_basis_set%set_radius(:) = set_radius(:)
    1800       31992 :       IF (PRESENT(lmax)) gto_basis_set%lmax(:) = lmax(:)
    1801       31992 :       IF (PRESENT(lmin)) gto_basis_set%lmin(:) = lmin(:)
    1802       31992 :       IF (PRESENT(lx)) gto_basis_set%lx(:) = lx(:)
    1803       31992 :       IF (PRESENT(ly)) gto_basis_set%ly(:) = ly(:)
    1804       31992 :       IF (PRESENT(lz)) gto_basis_set%lz(:) = lz(:)
    1805       31992 :       IF (PRESENT(m)) gto_basis_set%m(:) = m(:)
    1806       31992 :       IF (PRESENT(ncgf_set)) gto_basis_set%ncgf_set(:) = ncgf_set(:)
    1807       31992 :       IF (PRESENT(npgf)) gto_basis_set%npgf(:) = npgf(:)
    1808       31992 :       IF (PRESENT(nsgf_set)) gto_basis_set%nsgf_set(:) = nsgf_set(:)
    1809       31992 :       IF (PRESENT(nshell)) gto_basis_set%nshell(:) = nshell(:)
    1810       31992 :       IF (PRESENT(cphi)) gto_basis_set%cphi(:, :) = cphi(:, :)
    1811      434446 :       IF (PRESENT(pgf_radius)) gto_basis_set%pgf_radius(:, :) = pgf_radius(:, :)
    1812       31992 :       IF (PRESENT(sphi)) gto_basis_set%sphi(:, :) = sphi(:, :)
    1813       31992 :       IF (PRESENT(scon)) gto_basis_set%scon(:, :) = scon(:, :)
    1814       31992 :       IF (PRESENT(zet)) gto_basis_set%zet(:, :) = zet(:, :)
    1815       31992 :       IF (PRESENT(first_cgf)) gto_basis_set%first_cgf(:, :) = first_cgf(:, :)
    1816       31992 :       IF (PRESENT(first_sgf)) gto_basis_set%first_sgf(:, :) = first_sgf(:, :)
    1817       31992 :       IF (PRESENT(l)) l(:, :) = gto_basis_set%l(:, :)
    1818       31992 :       IF (PRESENT(last_cgf)) gto_basis_set%last_cgf(:, :) = last_cgf(:, :)
    1819       31992 :       IF (PRESENT(last_sgf)) gto_basis_set%last_sgf(:, :) = last_sgf(:, :)
    1820       31992 :       IF (PRESENT(n)) gto_basis_set%n(:, :) = n(:, :)
    1821       31992 :       IF (PRESENT(gcc)) gto_basis_set%gcc(:, :, :) = gcc(:, :, :)
    1822             : 
    1823       31992 :    END SUBROUTINE set_gto_basis_set
    1824             : 
    1825             : ! **************************************************************************************************
    1826             : !> \brief Write a Gaussian-type orbital (GTO) basis set data set to the output unit.
    1827             : !> \param gto_basis_set ...
    1828             : !> \param output_unit ...
    1829             : !> \param header ...
    1830             : !> \author MK
    1831             : ! **************************************************************************************************
    1832         139 :    SUBROUTINE write_gto_basis_set(gto_basis_set, output_unit, header)
    1833             : 
    1834             :       TYPE(gto_basis_set_type), INTENT(IN)               :: gto_basis_set
    1835             :       INTEGER, INTENT(in)                                :: output_unit
    1836             :       CHARACTER(len=*), OPTIONAL                         :: header
    1837             : 
    1838             :       INTEGER                                            :: ipgf, iset, ishell
    1839             : 
    1840         139 :       IF (output_unit > 0) THEN
    1841             : 
    1842         139 :          IF (PRESENT(header)) THEN
    1843             :             WRITE (UNIT=output_unit, FMT="(/,T6,A,T41,A40)") &
    1844         139 :                TRIM(header), TRIM(gto_basis_set%name)
    1845             :          END IF
    1846             : 
    1847             :          WRITE (UNIT=output_unit, FMT="(/,(T8,A,T71,I10))") &
    1848         139 :             "Number of orbital shell sets:            ", &
    1849         139 :             gto_basis_set%nset, &
    1850         139 :             "Number of orbital shells:                ", &
    1851         489 :             SUM(gto_basis_set%nshell(:)), &
    1852         139 :             "Number of primitive Cartesian functions: ", &
    1853         489 :             SUM(gto_basis_set%npgf(:)), &
    1854         139 :             "Number of Cartesian basis functions:     ", &
    1855         139 :             gto_basis_set%ncgf, &
    1856         139 :             "Number of spherical basis functions:     ", &
    1857         139 :             gto_basis_set%nsgf, &
    1858         139 :             "Norm type:                               ", &
    1859         278 :             gto_basis_set%norm_type
    1860             : 
    1861             :          WRITE (UNIT=output_unit, FMT="(/,T6,A,T41,A40,/,/,T25,A)") &
    1862         139 :             "GTO basis set information for", TRIM(gto_basis_set%name), &
    1863         278 :             "Set   Shell     n   l            Exponent    Coefficient"
    1864             : 
    1865         489 :          DO iset = 1, gto_basis_set%nset
    1866         350 :             WRITE (UNIT=output_unit, FMT="(A)") ""
    1867         976 :             DO ishell = 1, gto_basis_set%nshell(iset)
    1868             :                WRITE (UNIT=output_unit, &
    1869             :                       FMT="(T25,I3,4X,I4,4X,I2,2X,I2,(T51,2F15.6))") &
    1870         487 :                   iset, ishell, &
    1871         487 :                   gto_basis_set%n(ishell, iset), &
    1872         487 :                   gto_basis_set%l(ishell, iset), &
    1873        1839 :                   (gto_basis_set%zet(ipgf, iset), &
    1874        2326 :                    gto_basis_set%gcc(ipgf, ishell, iset), &
    1875        3650 :                    ipgf=1, gto_basis_set%npgf(iset))
    1876             :             END DO
    1877             :          END DO
    1878             : 
    1879             :       END IF
    1880             : 
    1881         139 :    END SUBROUTINE write_gto_basis_set
    1882             : 
    1883             : ! **************************************************************************************************
    1884             : 
    1885             : ! **************************************************************************************************
    1886             : !> \brief Write a Gaussian-type orbital (GTO) basis set data set to the output unit.
    1887             : !> \param orb_basis_set ...
    1888             : !> \param output_unit ...
    1889             : !> \param header ...
    1890             : !> \author MK
    1891             : ! **************************************************************************************************
    1892        4376 :    SUBROUTINE write_orb_basis_set(orb_basis_set, output_unit, header)
    1893             : 
    1894             :       TYPE(gto_basis_set_type), INTENT(IN)               :: orb_basis_set
    1895             :       INTEGER, INTENT(in)                                :: output_unit
    1896             :       CHARACTER(len=*), OPTIONAL                         :: header
    1897             : 
    1898             :       INTEGER                                            :: icgf, ico, ipgf, iset, ishell
    1899             : 
    1900        4376 :       IF (output_unit > 0) THEN
    1901        4376 :          IF (PRESENT(header)) THEN
    1902             :             WRITE (UNIT=output_unit, FMT="(/,T6,A,T41,A40)") &
    1903        4376 :                TRIM(header), TRIM(orb_basis_set%name)
    1904             :          END IF
    1905             : 
    1906             :          WRITE (UNIT=output_unit, FMT="(/,(T8,A,T71,I10))") &
    1907        4376 :             "Number of orbital shell sets:            ", &
    1908        4376 :             orb_basis_set%nset, &
    1909        4376 :             "Number of orbital shells:                ", &
    1910       16711 :             SUM(orb_basis_set%nshell(:)), &
    1911        4376 :             "Number of primitive Cartesian functions: ", &
    1912       16711 :             SUM(orb_basis_set%npgf(:)), &
    1913        4376 :             "Number of Cartesian basis functions:     ", &
    1914        4376 :             orb_basis_set%ncgf, &
    1915        4376 :             "Number of spherical basis functions:     ", &
    1916        4376 :             orb_basis_set%nsgf, &
    1917        4376 :             "Norm type:                               ", &
    1918        8752 :             orb_basis_set%norm_type
    1919             : 
    1920             :          WRITE (UNIT=output_unit, FMT="(/,T8,A,/,/,T25,A)") &
    1921        4376 :             "Normalised Cartesian orbitals:", &
    1922        8752 :             "Set   Shell   Orbital            Exponent    Coefficient"
    1923             : 
    1924        4376 :          icgf = 0
    1925             : 
    1926       16711 :          DO iset = 1, orb_basis_set%nset
    1927       35244 :             DO ishell = 1, orb_basis_set%nshell(iset)
    1928       18533 :                WRITE (UNIT=output_unit, FMT="(A)") ""
    1929       81190 :                DO ico = 1, nco(orb_basis_set%l(ishell, iset))
    1930       50322 :                   icgf = icgf + 1
    1931             :                   WRITE (UNIT=output_unit, &
    1932             :                          FMT="(T25,I3,4X,I4,3X,A12,(T51,2F15.6))") &
    1933       50322 :                      iset, ishell, orb_basis_set%cgf_symbol(icgf), &
    1934      157916 :                      (orb_basis_set%zet(ipgf, iset), &
    1935             :                       orb_basis_set%norm_cgf(icgf)* &
    1936      208238 :                       orb_basis_set%gcc(ipgf, ishell, iset), &
    1937      327415 :                       ipgf=1, orb_basis_set%npgf(iset))
    1938             :                END DO
    1939             :             END DO
    1940             :          END DO
    1941             :       END IF
    1942             : 
    1943        4376 :    END SUBROUTINE write_orb_basis_set
    1944             : 
    1945             : ! **************************************************************************************************
    1946             : !> \brief ...
    1947             : !> \param sto_basis_set ...
    1948             : ! **************************************************************************************************
    1949        3048 :    SUBROUTINE allocate_sto_basis_set(sto_basis_set)
    1950             : 
    1951             :       TYPE(sto_basis_set_type), POINTER                  :: sto_basis_set
    1952             : 
    1953        3048 :       CALL deallocate_sto_basis_set(sto_basis_set)
    1954             : 
    1955        3048 :       ALLOCATE (sto_basis_set)
    1956             : 
    1957        3048 :    END SUBROUTINE allocate_sto_basis_set
    1958             : 
    1959             : ! **************************************************************************************************
    1960             : !> \brief ...
    1961             : !> \param sto_basis_set ...
    1962             : ! **************************************************************************************************
    1963        7820 :    SUBROUTINE deallocate_sto_basis_set(sto_basis_set)
    1964             : 
    1965             :       TYPE(sto_basis_set_type), POINTER                  :: sto_basis_set
    1966             : 
    1967             : ! -------------------------------------------------------------------------
    1968             : 
    1969        7820 :       IF (ASSOCIATED(sto_basis_set)) THEN
    1970        3048 :          IF (ASSOCIATED(sto_basis_set%symbol)) THEN
    1971        2908 :             DEALLOCATE (sto_basis_set%symbol)
    1972             :          END IF
    1973        3048 :          IF (ASSOCIATED(sto_basis_set%nq)) THEN
    1974        3048 :             DEALLOCATE (sto_basis_set%nq)
    1975             :          END IF
    1976        3048 :          IF (ASSOCIATED(sto_basis_set%lq)) THEN
    1977        3048 :             DEALLOCATE (sto_basis_set%lq)
    1978             :          END IF
    1979        3048 :          IF (ASSOCIATED(sto_basis_set%zet)) THEN
    1980        3048 :             DEALLOCATE (sto_basis_set%zet)
    1981             :          END IF
    1982             : 
    1983        3048 :          DEALLOCATE (sto_basis_set)
    1984             :       END IF
    1985        7820 :    END SUBROUTINE deallocate_sto_basis_set
    1986             : 
    1987             : ! **************************************************************************************************
    1988             : !> \brief ...
    1989             : !> \param sto_basis_set ...
    1990             : !> \param name ...
    1991             : !> \param nshell ...
    1992             : !> \param symbol ...
    1993             : !> \param nq ...
    1994             : !> \param lq ...
    1995             : !> \param zet ...
    1996             : !> \param maxlq ...
    1997             : !> \param numsto ...
    1998             : ! **************************************************************************************************
    1999        3048 :    SUBROUTINE get_sto_basis_set(sto_basis_set, name, nshell, symbol, nq, lq, zet, maxlq, numsto)
    2000             : 
    2001             :       TYPE(sto_basis_set_type), INTENT(IN)               :: sto_basis_set
    2002             :       CHARACTER(LEN=default_string_length), &
    2003             :          INTENT(OUT), OPTIONAL                           :: name
    2004             :       INTEGER, INTENT(OUT), OPTIONAL                     :: nshell
    2005             :       CHARACTER(LEN=6), DIMENSION(:), OPTIONAL, POINTER  :: symbol
    2006             :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: nq, lq
    2007             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: zet
    2008             :       INTEGER, INTENT(OUT), OPTIONAL                     :: maxlq, numsto
    2009             : 
    2010             :       INTEGER                                            :: iset
    2011             : 
    2012        3048 :       IF (PRESENT(name)) name = sto_basis_set%name
    2013        3048 :       IF (PRESENT(nshell)) nshell = sto_basis_set%nshell
    2014        3048 :       IF (PRESENT(symbol)) symbol => sto_basis_set%symbol
    2015        3048 :       IF (PRESENT(nq)) nq => sto_basis_set%nq
    2016        3048 :       IF (PRESENT(lq)) lq => sto_basis_set%lq
    2017        3048 :       IF (PRESENT(zet)) zet => sto_basis_set%zet
    2018        3048 :       IF (PRESENT(maxlq)) THEN
    2019           0 :          maxlq = MAXVAL(sto_basis_set%lq(1:sto_basis_set%nshell))
    2020             :       END IF
    2021        3048 :       IF (PRESENT(numsto)) THEN
    2022           0 :          numsto = 0
    2023           0 :          DO iset = 1, sto_basis_set%nshell
    2024           0 :             numsto = numsto + 2*sto_basis_set%lq(iset) + 1
    2025             :          END DO
    2026             :       END IF
    2027             : 
    2028        3048 :    END SUBROUTINE get_sto_basis_set
    2029             : 
    2030             : ! **************************************************************************************************
    2031             : !> \brief ...
    2032             : !> \param sto_basis_set ...
    2033             : !> \param name ...
    2034             : !> \param nshell ...
    2035             : !> \param symbol ...
    2036             : !> \param nq ...
    2037             : !> \param lq ...
    2038             : !> \param zet ...
    2039             : ! **************************************************************************************************
    2040        3044 :    SUBROUTINE set_sto_basis_set(sto_basis_set, name, nshell, symbol, nq, lq, zet)
    2041             : 
    2042             :       TYPE(sto_basis_set_type), INTENT(INOUT)            :: sto_basis_set
    2043             :       CHARACTER(LEN=default_string_length), INTENT(IN), &
    2044             :          OPTIONAL                                        :: name
    2045             :       INTEGER, INTENT(IN), OPTIONAL                      :: nshell
    2046             :       CHARACTER(LEN=6), DIMENSION(:), OPTIONAL, POINTER  :: symbol
    2047             :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: nq, lq
    2048             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: zet
    2049             : 
    2050             :       INTEGER                                            :: ns
    2051             : 
    2052        3044 :       IF (PRESENT(name)) sto_basis_set%name = name
    2053        3044 :       IF (PRESENT(nshell)) sto_basis_set%nshell = nshell
    2054        3044 :       IF (PRESENT(symbol)) THEN
    2055        2904 :          ns = SIZE(symbol)
    2056        2904 :          IF (ASSOCIATED(sto_basis_set%symbol)) DEALLOCATE (sto_basis_set%symbol)
    2057        8712 :          ALLOCATE (sto_basis_set%symbol(1:ns))
    2058       13240 :          sto_basis_set%symbol(:) = symbol(:)
    2059             :       END IF
    2060        3044 :       IF (PRESENT(nq)) THEN
    2061        3044 :          ns = SIZE(nq)
    2062        3044 :          CALL reallocate(sto_basis_set%nq, 1, ns)
    2063       13520 :          sto_basis_set%nq = nq(:)
    2064             :       END IF
    2065        3044 :       IF (PRESENT(lq)) THEN
    2066        3044 :          ns = SIZE(lq)
    2067        3044 :          CALL reallocate(sto_basis_set%lq, 1, ns)
    2068       13520 :          sto_basis_set%lq = lq(:)
    2069             :       END IF
    2070        3044 :       IF (PRESENT(zet)) THEN
    2071        3044 :          ns = SIZE(zet)
    2072        3044 :          CALL reallocate(sto_basis_set%zet, 1, ns)
    2073       13520 :          sto_basis_set%zet = zet(:)
    2074             :       END IF
    2075             : 
    2076        3044 :    END SUBROUTINE set_sto_basis_set
    2077             : 
    2078             : ! **************************************************************************************************
    2079             : !> \brief ...
    2080             : !> \param element_symbol ...
    2081             : !> \param basis_set_name ...
    2082             : !> \param sto_basis_set ...
    2083             : !> \param para_env ...
    2084             : !> \param dft_section ...
    2085             : ! **************************************************************************************************
    2086           4 :    SUBROUTINE read_sto_basis_set(element_symbol, basis_set_name, sto_basis_set, para_env, dft_section)
    2087             : 
    2088             :       ! Read a Slater-type orbital (STO) basis set from the database file.
    2089             : 
    2090             :       CHARACTER(LEN=*), INTENT(IN)                       :: element_symbol, basis_set_name
    2091             :       TYPE(sto_basis_set_type), INTENT(INOUT)            :: sto_basis_set
    2092             :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2093             :       TYPE(section_vals_type), POINTER                   :: dft_section
    2094             : 
    2095             :       CHARACTER(LEN=10)                                  :: nlsym
    2096             :       CHARACTER(LEN=2)                                   :: lsym
    2097             :       CHARACTER(LEN=240)                                 :: line
    2098             :       CHARACTER(LEN=242)                                 :: line2
    2099             :       CHARACTER(len=default_path_length)                 :: basis_set_file_name, tmp
    2100             :       CHARACTER(LEN=default_path_length), DIMENSION(:), &
    2101           4 :          POINTER                                         :: cbasis
    2102           4 :       CHARACTER(LEN=LEN(basis_set_name))                 :: bsname
    2103           4 :       CHARACTER(LEN=LEN(basis_set_name)+2)               :: bsname2
    2104           4 :       CHARACTER(LEN=LEN(element_symbol))                 :: symbol
    2105           4 :       CHARACTER(LEN=LEN(element_symbol)+2)               :: symbol2
    2106             :       INTEGER                                            :: ibasis, irep, iset, nbasis, nq, nset, &
    2107             :                                                             strlen1, strlen2
    2108             :       LOGICAL                                            :: basis_found, found, match
    2109             :       REAL(KIND=dp)                                      :: zet
    2110             :       TYPE(cp_parser_type)                               :: parser
    2111             : 
    2112           4 :       line = ""
    2113           4 :       line2 = ""
    2114           4 :       symbol = ""
    2115           4 :       symbol2 = ""
    2116           4 :       bsname = ""
    2117           4 :       bsname2 = ""
    2118             : 
    2119           4 :       nbasis = 1
    2120             : 
    2121           4 :       sto_basis_set%name = basis_set_name
    2122             :       CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", &
    2123           4 :                                 n_rep_val=nbasis)
    2124          12 :       ALLOCATE (cbasis(nbasis))
    2125           8 :       DO ibasis = 1, nbasis
    2126             :          CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", &
    2127           4 :                                    i_rep_val=ibasis, c_val=cbasis(ibasis))
    2128           4 :          basis_set_file_name = cbasis(ibasis)
    2129           4 :          tmp = basis_set_file_name
    2130           8 :          CALL uppercase(tmp)
    2131             :       END DO
    2132             : 
    2133             :       ! Search for the requested basis set in the basis set file
    2134             :       ! until the basis set is found or the end of file is reached
    2135             : 
    2136           4 :       basis_found = .FALSE.
    2137           8 :       basis_loop: DO ibasis = 1, nbasis
    2138           4 :          IF (basis_found) EXIT basis_loop
    2139           4 :          basis_set_file_name = cbasis(ibasis)
    2140           4 :          CALL parser_create(parser, basis_set_file_name, para_env=para_env)
    2141             : 
    2142           4 :          bsname = basis_set_name
    2143           4 :          symbol = element_symbol
    2144           4 :          irep = 0
    2145             : 
    2146           4 :          tmp = basis_set_name
    2147           4 :          CALL uppercase(tmp)
    2148             : 
    2149           4 :          IF (tmp .NE. "NONE") THEN
    2150             :             search_loop: DO
    2151             : 
    2152           6 :                CALL parser_search_string(parser, TRIM(bsname), .TRUE., found, line)
    2153           6 :                IF (found) THEN
    2154           6 :                   CALL uppercase(symbol)
    2155           6 :                   CALL uppercase(bsname)
    2156             : 
    2157           6 :                   match = .FALSE.
    2158           6 :                   CALL uppercase(line)
    2159             :                   ! Check both the element symbol and the basis set name
    2160           6 :                   line2 = " "//line//" "
    2161           6 :                   symbol2 = " "//TRIM(symbol)//" "
    2162           6 :                   bsname2 = " "//TRIM(bsname)//" "
    2163           6 :                   strlen1 = LEN_TRIM(symbol2) + 1
    2164           6 :                   strlen2 = LEN_TRIM(bsname2) + 1
    2165             : 
    2166           6 :                   IF ((INDEX(line2, symbol2(:strlen1)) > 0) .AND. &
    2167             :                       (INDEX(line2, bsname2(:strlen2)) > 0)) match = .TRUE.
    2168             :                   IF (match) THEN
    2169             :                      ! Read the basis set information
    2170           4 :                      CALL parser_get_object(parser, nset, newline=.TRUE.)
    2171           4 :                      sto_basis_set%nshell = nset
    2172             : 
    2173           4 :                      CALL reallocate(sto_basis_set%nq, 1, nset)
    2174           4 :                      CALL reallocate(sto_basis_set%lq, 1, nset)
    2175           4 :                      CALL reallocate(sto_basis_set%zet, 1, nset)
    2176          12 :                      ALLOCATE (sto_basis_set%symbol(nset))
    2177             : 
    2178          20 :                      DO iset = 1, nset
    2179          16 :                         CALL parser_get_object(parser, nq, newline=.TRUE.)
    2180          16 :                         CALL parser_get_object(parser, lsym)
    2181          16 :                         CALL parser_get_object(parser, zet)
    2182          16 :                         sto_basis_set%nq(iset) = nq
    2183          16 :                         sto_basis_set%zet(iset) = zet
    2184          16 :                         WRITE (nlsym, "(I2,A)") nq, TRIM(lsym)
    2185          16 :                         sto_basis_set%symbol(iset) = TRIM(nlsym)
    2186          20 :                         SELECT CASE (TRIM(lsym))
    2187             :                         CASE ("S", "s")
    2188          12 :                            sto_basis_set%lq(iset) = 0
    2189             :                         CASE ("P", "p")
    2190           4 :                            sto_basis_set%lq(iset) = 1
    2191             :                         CASE ("D", "d")
    2192           0 :                            sto_basis_set%lq(iset) = 2
    2193             :                         CASE ("F", "f")
    2194           0 :                            sto_basis_set%lq(iset) = 3
    2195             :                         CASE ("G", "g")
    2196           0 :                            sto_basis_set%lq(iset) = 4
    2197             :                         CASE ("H", "h")
    2198           0 :                            sto_basis_set%lq(iset) = 5
    2199             :                         CASE ("I", "i", "J", "j")
    2200           0 :                            sto_basis_set%lq(iset) = 6
    2201             :                         CASE ("K", "k")
    2202           0 :                            sto_basis_set%lq(iset) = 7
    2203             :                         CASE ("L", "l")
    2204           0 :                            sto_basis_set%lq(iset) = 8
    2205             :                         CASE ("M", "m")
    2206           0 :                            sto_basis_set%lq(iset) = 9
    2207             :                         CASE DEFAULT
    2208             :                            CALL cp_abort(__LOCATION__, &
    2209             :                                          "The requested basis set <"//TRIM(bsname)// &
    2210          16 :                                          "> for element <"//TRIM(symbol)//"> has an invalid component: ")
    2211             :                         END SELECT
    2212             :                      END DO
    2213             : 
    2214             :                      basis_found = .TRUE.
    2215             :                      EXIT search_loop
    2216             :                   END IF
    2217             :                ELSE
    2218             :                   EXIT search_loop
    2219             :                END IF
    2220             :             END DO search_loop
    2221             :          ELSE
    2222           4 :             match = .FALSE.
    2223             :          END IF
    2224             : 
    2225           8 :          CALL parser_release(parser)
    2226             : 
    2227             :       END DO basis_loop
    2228             : 
    2229           4 :       IF (tmp .NE. "NONE") THEN
    2230           4 :          IF (.NOT. basis_found) THEN
    2231           0 :             basis_set_file_name = ""
    2232           0 :             DO ibasis = 1, nbasis
    2233           0 :                basis_set_file_name = TRIM(basis_set_file_name)//"<"//TRIM(cbasis(ibasis))//"> "
    2234             :             END DO
    2235             :             CALL cp_abort(__LOCATION__, &
    2236             :                           "The requested basis set <"//TRIM(bsname)// &
    2237             :                           "> for element <"//TRIM(symbol)//"> was not "// &
    2238             :                           "found in the basis set files "// &
    2239           0 :                           TRIM(basis_set_file_name))
    2240             :          END IF
    2241             :       END IF
    2242           4 :       DEALLOCATE (cbasis)
    2243             : 
    2244          16 :    END SUBROUTINE read_sto_basis_set
    2245             : 
    2246             : ! **************************************************************************************************
    2247             : !> \brief ...
    2248             : !> \param sto_basis_set ...
    2249             : !> \param gto_basis_set ...
    2250             : !> \param ngauss ...
    2251             : !> \param ortho ...
    2252             : ! **************************************************************************************************
    2253        3048 :    SUBROUTINE create_gto_from_sto_basis(sto_basis_set, gto_basis_set, ngauss, ortho)
    2254             : 
    2255             :       TYPE(sto_basis_set_type), INTENT(IN)               :: sto_basis_set
    2256             :       TYPE(gto_basis_set_type), POINTER                  :: gto_basis_set
    2257             :       INTEGER, INTENT(IN), OPTIONAL                      :: ngauss
    2258             :       LOGICAL, INTENT(IN), OPTIONAL                      :: ortho
    2259             : 
    2260             :       INTEGER, PARAMETER                                 :: maxng = 6
    2261             : 
    2262             :       CHARACTER(LEN=default_string_length)               :: name, sng
    2263             :       INTEGER                                            :: i1, i2, ico, ipgf, iset, jset, l, &
    2264             :                                                             lshell, m, maxco, maxl, ncgf, ng, ngs, &
    2265             :                                                             np, nset, nsgf, nshell
    2266             :       INTEGER, DIMENSION(0:10)                           :: mxf
    2267        3048 :       INTEGER, DIMENSION(:), POINTER                     :: lq, nq
    2268             :       LOGICAL                                            :: do_ortho
    2269        3048 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: gal, zal, zll
    2270        3048 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: zet
    2271             :       REAL(KIND=dp), DIMENSION(maxng)                    :: gcc, zetg
    2272             : 
    2273        3048 :       ng = 6
    2274        3048 :       IF (PRESENT(ngauss)) ng = ngauss
    2275        3048 :       IF (ng > maxng) CPABORT("Too many Gaussian primitives requested")
    2276        3048 :       do_ortho = .FALSE.
    2277        3048 :       IF (PRESENT(ortho)) do_ortho = ortho
    2278             : 
    2279        3048 :       CALL allocate_gto_basis_set(gto_basis_set)
    2280             : 
    2281             :       CALL get_sto_basis_set(sto_basis_set, name=name, nshell=nshell, nq=nq, &
    2282        3048 :                              lq=lq, zet=zet)
    2283             : 
    2284       13540 :       maxl = MAXVAL(lq)
    2285        3048 :       CALL init_orbital_pointers(maxl)
    2286             : 
    2287        3048 :       CALL integer_to_string(ng, sng)
    2288        3048 :       gto_basis_set%name = TRIM(name)//"_STO-"//TRIM(sng)//"G"
    2289             : 
    2290        3048 :       nset = nshell
    2291        3048 :       gto_basis_set%nset = nset
    2292        3048 :       CALL reallocate(gto_basis_set%lmax, 1, nset)
    2293        3048 :       CALL reallocate(gto_basis_set%lmin, 1, nset)
    2294        3048 :       CALL reallocate(gto_basis_set%npgf, 1, nset)
    2295        3048 :       CALL reallocate(gto_basis_set%nshell, 1, nset)
    2296        3048 :       CALL reallocate(gto_basis_set%n, 1, 1, 1, nset)
    2297        3048 :       CALL reallocate(gto_basis_set%l, 1, 1, 1, nset)
    2298        3048 :       CALL reallocate(gto_basis_set%zet, 1, ng, 1, nset)
    2299        3048 :       CALL reallocate(gto_basis_set%gcc, 1, ng, 1, 1, 1, nset)
    2300             : 
    2301        9330 :       DO iset = 1, nset
    2302        6282 :          CALL get_sto_ng(zet(iset), ng, nq(iset), lq(iset), zetg, gcc)
    2303        6282 :          gto_basis_set%lmax(iset) = lq(iset)
    2304        6282 :          gto_basis_set%lmin(iset) = lq(iset)
    2305        6282 :          gto_basis_set%npgf(iset) = ng
    2306        6282 :          gto_basis_set%nshell(iset) = 1
    2307        6282 :          gto_basis_set%n(1, iset) = lq(iset) + 1
    2308        6282 :          gto_basis_set%l(1, iset) = lq(iset)
    2309       45300 :          DO ipgf = 1, ng
    2310       35970 :             gto_basis_set%gcc(ipgf, 1, iset) = gcc(ipgf)
    2311       42252 :             gto_basis_set%zet(ipgf, iset) = zetg(ipgf)
    2312             :          END DO
    2313             :       END DO
    2314             : 
    2315        3048 :       IF (do_ortho) THEN
    2316         590 :          mxf = 0
    2317        1844 :          DO iset = 1, nset
    2318        1254 :             l = gto_basis_set%l(1, iset)
    2319        1844 :             mxf(l) = mxf(l) + 1
    2320             :          END DO
    2321        7080 :          m = MAXVAL(mxf)
    2322         590 :          IF (m > 1) THEN
    2323        2260 :             ALLOCATE (gal(ng, nset), zal(ng, nset), zll(m*ng, 0:maxl))
    2324         678 :             DO iset = 1, nset
    2325        2260 :                zal(1:ng, iset) = gto_basis_set%zet(1:ng, iset)
    2326        2486 :                gal(1:ng, iset) = gto_basis_set%gcc(1:ng, 1, iset)
    2327             :             END DO
    2328         226 :             CALL reallocate(gto_basis_set%zet, 1, m*ng, 1, nset)
    2329         226 :             CALL reallocate(gto_basis_set%gcc, 1, m*ng, 1, 1, 1, nset)
    2330         678 :             DO iset = 1, nset
    2331         452 :                l = gto_basis_set%l(1, iset)
    2332         678 :                gto_basis_set%npgf(iset) = ng*mxf(l)
    2333             :             END DO
    2334        4294 :             gto_basis_set%zet = 0.0_dp
    2335        4746 :             gto_basis_set%gcc = 0.0_dp
    2336        2260 :             zll = 0.0_dp
    2337         226 :             mxf = 0
    2338         678 :             DO iset = 1, nset
    2339         452 :                l = gto_basis_set%l(1, iset)
    2340         452 :                mxf(l) = mxf(l) + 1
    2341         452 :                i1 = mxf(l)*ng - ng + 1
    2342         452 :                i2 = mxf(l)*ng
    2343        2260 :                zll(i1:i2, l) = zal(1:ng, iset)
    2344        2486 :                gto_basis_set%gcc(i1:i2, 1, iset) = gal(1:ng, iset)
    2345             :             END DO
    2346         678 :             DO iset = 1, nset
    2347         452 :                l = gto_basis_set%l(1, iset)
    2348        4294 :                gto_basis_set%zet(:, iset) = zll(:, l)
    2349             :             END DO
    2350         678 :             DO iset = 1, nset
    2351         452 :                l = gto_basis_set%l(1, iset)
    2352         904 :                DO jset = 1, iset - 1
    2353         678 :                   IF (gto_basis_set%l(1, iset) == l) THEN
    2354         226 :                      m = mxf(l)*ng
    2355             :                      CALL orthofun(gto_basis_set%zet(1:m, iset), gto_basis_set%gcc(1:m, 1, iset), &
    2356         226 :                                    gto_basis_set%gcc(1:m, 1, jset), l)
    2357             :                   END IF
    2358             :                END DO
    2359             :             END DO
    2360         226 :             DEALLOCATE (gal, zal, zll)
    2361             :          END IF
    2362             :       END IF
    2363             : 
    2364        9330 :       ngs = MAXVAL(gto_basis_set%npgf(1:nset))
    2365        3048 :       CALL reallocate(gto_basis_set%set_radius, 1, nset)
    2366        3048 :       CALL reallocate(gto_basis_set%pgf_radius, 1, ngs, 1, nset)
    2367        3048 :       CALL reallocate(gto_basis_set%first_cgf, 1, 1, 1, nset)
    2368        3048 :       CALL reallocate(gto_basis_set%first_sgf, 1, 1, 1, nset)
    2369        3048 :       CALL reallocate(gto_basis_set%last_cgf, 1, 1, 1, nset)
    2370        3048 :       CALL reallocate(gto_basis_set%last_sgf, 1, 1, 1, nset)
    2371        3048 :       CALL reallocate(gto_basis_set%ncgf_set, 1, nset)
    2372        3048 :       CALL reallocate(gto_basis_set%nsgf_set, 1, nset)
    2373             : 
    2374        3048 :       maxco = 0
    2375        3048 :       ncgf = 0
    2376        3048 :       nsgf = 0
    2377             : 
    2378        9330 :       DO iset = 1, nset
    2379        6282 :          gto_basis_set%ncgf_set(iset) = 0
    2380        6282 :          gto_basis_set%nsgf_set(iset) = 0
    2381        6282 :          lshell = gto_basis_set%l(1, iset)
    2382        6282 :          gto_basis_set%first_cgf(1, iset) = ncgf + 1
    2383        6282 :          ncgf = ncgf + nco(lshell)
    2384        6282 :          gto_basis_set%last_cgf(1, iset) = ncgf
    2385             :          gto_basis_set%ncgf_set(iset) = &
    2386        6282 :             gto_basis_set%ncgf_set(iset) + nco(lshell)
    2387        6282 :          gto_basis_set%first_sgf(1, iset) = nsgf + 1
    2388        6282 :          nsgf = nsgf + nso(lshell)
    2389        6282 :          gto_basis_set%last_sgf(1, iset) = nsgf
    2390             :          gto_basis_set%nsgf_set(iset) = &
    2391        6282 :             gto_basis_set%nsgf_set(iset) + nso(lshell)
    2392        6282 :          ngs = gto_basis_set%npgf(iset)
    2393        9330 :          maxco = MAX(maxco, ngs*ncoset(lshell))
    2394             :       END DO
    2395             : 
    2396        3048 :       gto_basis_set%ncgf = ncgf
    2397        3048 :       gto_basis_set%nsgf = nsgf
    2398             : 
    2399        3048 :       CALL reallocate(gto_basis_set%cphi, 1, maxco, 1, ncgf)
    2400        3048 :       CALL reallocate(gto_basis_set%sphi, 1, maxco, 1, nsgf)
    2401        3048 :       CALL reallocate(gto_basis_set%scon, 1, maxco, 1, nsgf)
    2402        3048 :       CALL reallocate(gto_basis_set%lx, 1, ncgf)
    2403        3048 :       CALL reallocate(gto_basis_set%ly, 1, ncgf)
    2404        3048 :       CALL reallocate(gto_basis_set%lz, 1, ncgf)
    2405        3048 :       CALL reallocate(gto_basis_set%m, 1, nsgf)
    2406        3048 :       CALL reallocate(gto_basis_set%norm_cgf, 1, ncgf)
    2407        9144 :       ALLOCATE (gto_basis_set%cgf_symbol(ncgf))
    2408        9144 :       ALLOCATE (gto_basis_set%sgf_symbol(nsgf))
    2409             : 
    2410        3048 :       ncgf = 0
    2411        3048 :       nsgf = 0
    2412             : 
    2413        9330 :       DO iset = 1, nset
    2414        6282 :          lshell = gto_basis_set%l(1, iset)
    2415        6282 :          np = lshell + 1
    2416       21004 :          DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
    2417       14722 :             ncgf = ncgf + 1
    2418       14722 :             gto_basis_set%lx(ncgf) = indco(1, ico)
    2419       14722 :             gto_basis_set%ly(ncgf) = indco(2, ico)
    2420       14722 :             gto_basis_set%lz(ncgf) = indco(3, ico)
    2421             :             gto_basis_set%cgf_symbol(ncgf) = &
    2422             :                cgf_symbol(np, (/gto_basis_set%lx(ncgf), &
    2423             :                                 gto_basis_set%ly(ncgf), &
    2424       65170 :                                 gto_basis_set%lz(ncgf)/))
    2425             :          END DO
    2426       23240 :          DO m = -lshell, lshell
    2427       13910 :             nsgf = nsgf + 1
    2428       13910 :             gto_basis_set%m(nsgf) = m
    2429       20192 :             gto_basis_set%sgf_symbol(nsgf) = sgf_symbol(np, lshell, m)
    2430             :          END DO
    2431             :       END DO
    2432             : 
    2433        3048 :       gto_basis_set%norm_type = -1
    2434             : 
    2435        6096 :    END SUBROUTINE create_gto_from_sto_basis
    2436             : 
    2437             : ! **************************************************************************************************
    2438             : !> \brief ...
    2439             : !> \param zet ...
    2440             : !> \param co ...
    2441             : !> \param cr ...
    2442             : !> \param l ...
    2443             : ! **************************************************************************************************
    2444         226 :    SUBROUTINE orthofun(zet, co, cr, l)
    2445             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zet
    2446             :       REAL(KIND=dp), DIMENSION(:), INTENT(INOUT)         :: co, cr
    2447             :       INTEGER, INTENT(IN)                                :: l
    2448             : 
    2449             :       REAL(KIND=dp)                                      :: ss
    2450             : 
    2451         226 :       CALL aovlp(l, zet, cr, cr, ss)
    2452        2034 :       cr(:) = cr(:)/SQRT(ss)
    2453         226 :       CALL aovlp(l, zet, co, cr, ss)
    2454        2034 :       co(:) = co(:) - ss*cr(:)
    2455         226 :       CALL aovlp(l, zet, co, co, ss)
    2456        2034 :       co(:) = co(:)/SQRT(ss)
    2457             : 
    2458         226 :    END SUBROUTINE orthofun
    2459             : 
    2460             : ! **************************************************************************************************
    2461             : !> \brief ...
    2462             : !> \param l ...
    2463             : !> \param zet ...
    2464             : !> \param ca ...
    2465             : !> \param cb ...
    2466             : !> \param ss ...
    2467             : ! **************************************************************************************************
    2468         678 :    SUBROUTINE aovlp(l, zet, ca, cb, ss)
    2469             :       INTEGER, INTENT(IN)                                :: l
    2470             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN)            :: zet, ca, cb
    2471             :       REAL(KIND=dp), INTENT(OUT)                         :: ss
    2472             : 
    2473             :       INTEGER                                            :: i, j, m
    2474             :       REAL(KIND=dp)                                      :: ab, ai, aj, s00, sss
    2475             : 
    2476             : !
    2477             : ! use init_norm_cgf_orb
    2478             : !
    2479         678 :       m = SIZE(zet)
    2480         678 :       ss = 0.0_dp
    2481        6102 :       DO i = 1, m
    2482        5424 :          ai = (2.0_dp*zet(i)/pi)**0.75_dp
    2483       49494 :          DO j = 1, m
    2484       43392 :             aj = (2.0_dp*zet(j)/pi)**0.75_dp
    2485       43392 :             ab = 1._dp/(zet(i) + zet(j))
    2486       43392 :             s00 = ai*aj*(pi*ab)**1.50_dp
    2487       43392 :             IF (l == 0) THEN
    2488             :                sss = s00
    2489           0 :             ELSEIF (l == 1) THEN
    2490           0 :                sss = s00*ab*0.5_dp
    2491             :             ELSE
    2492           0 :                CPABORT("aovlp lvalue")
    2493             :             END IF
    2494       48816 :             ss = ss + sss*ca(i)*cb(j)
    2495             :          END DO
    2496             :       END DO
    2497             : 
    2498         678 :    END SUBROUTINE aovlp
    2499             : 
    2500             : ! **************************************************************************************************
    2501             : !> \brief ...
    2502             : !> \param z ...
    2503             : !> \param ne ...
    2504             : !> \param n ...
    2505             : !> \param l ...
    2506             : !> \return ...
    2507             : ! **************************************************************************************************
    2508       21284 :    PURE FUNCTION srules(z, ne, n, l)
    2509             :       ! Slater rules
    2510             :       INTEGER, INTENT(IN)                                :: z
    2511             :       INTEGER, DIMENSION(:, :), INTENT(IN)               :: ne
    2512             :       INTEGER, INTENT(IN)                                :: n, l
    2513             :       REAL(dp)                                           :: srules
    2514             : 
    2515             :       REAL(dp), DIMENSION(7), PARAMETER :: &
    2516             :          xns = (/1.0_dp, 2.0_dp, 3.0_dp, 3.7_dp, 4.0_dp, 4.2_dp, 4.4_dp/)
    2517             : 
    2518             :       INTEGER                                            :: i, l1, l2, m, m1, m2, nn
    2519             :       REAL(dp)                                           :: s
    2520             : 
    2521       21284 :       s = 0.0_dp
    2522             :       ! The complete shell
    2523       21284 :       l1 = MIN(l + 1, 4)
    2524       21284 :       nn = MIN(n, 7)
    2525             :       IF (l1 == 1) l2 = 2
    2526       21284 :       IF (l1 == 2) l2 = 1
    2527       21284 :       IF (l1 == 3) l2 = 4
    2528       21284 :       IF (l1 == 4) l2 = 3
    2529             :       ! Rule a) no contribution from shells further out
    2530             :       ! Rule b) 0.35 (1s 0.3) from each other electron in the same shell
    2531       21284 :       IF (n == 1) THEN
    2532        5962 :          m = ne(1, 1)
    2533        5962 :          s = s + 0.3_dp*REAL(m - 1, dp)
    2534             :       ELSE
    2535       15322 :          m = ne(l1, nn) + ne(l2, nn)
    2536       15322 :          s = s + 0.35_dp*REAL(m - 1, dp)
    2537             :       END IF
    2538             :       ! Rule c) if (s,p) shell 0.85 from each electron with n-1, and 1.0
    2539             :       !      from all electrons further in
    2540       21284 :       IF (l1 + l2 == 3) THEN
    2541       19738 :          IF (nn > 1) THEN
    2542       13776 :             m1 = ne(1, nn - 1) + ne(2, nn - 1) + ne(3, nn - 1) + ne(4, nn - 1)
    2543       13776 :             m2 = 0
    2544       21524 :             DO i = 1, nn - 2
    2545       21524 :                m2 = m2 + ne(1, i) + ne(2, i) + ne(3, i) + ne(4, I)
    2546             :             END DO
    2547       13776 :             s = s + 0.85_dp*REAL(m1, dp) + 1._dp*REAL(m2, dp)
    2548             :          END IF
    2549             :       ELSE
    2550             :          ! Rule d) if (d,f) shell 1.0 from each electron inside
    2551             :          m = 0
    2552        6274 :          DO i = 1, nn - 1
    2553        6274 :             m = m + ne(1, i) + ne(2, i) + ne(3, i) + ne(4, i)
    2554             :          END DO
    2555        1546 :          s = s + 1._dp*REAL(m, dp)
    2556             :       END IF
    2557             :       ! Slater exponent is (Z-S)/NS
    2558       21284 :       srules = (REAL(z, dp) - s)/xns(nn)
    2559       21284 :    END FUNCTION srules
    2560             : 
    2561             : ! **************************************************************************************************
    2562             : !> \brief sort basis sets w.r.t. radius
    2563             : !> \param basis_set ...
    2564             : !> \param sort_method ...
    2565             : ! **************************************************************************************************
    2566       11324 :    SUBROUTINE sort_gto_basis_set(basis_set, sort_method)
    2567             :       TYPE(gto_basis_set_type), INTENT(INOUT)            :: basis_set
    2568             :       INTEGER, INTENT(IN)                                :: sort_method
    2569             : 
    2570       11324 :       CHARACTER(LEN=12), DIMENSION(:), POINTER           :: cgf_symbol
    2571       11324 :       CHARACTER(LEN=6), DIMENSION(:), POINTER            :: sgf_symbol
    2572             :       INTEGER :: ic, ic_max, icgf, icgf_new, icgf_old, ico, is, is_max, iset, isgf, isgf_new, &
    2573             :          isgf_old, ishell, lshell, maxco, maxpgf, maxshell, mm, nc, ncgf, ns, nset
    2574       11324 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: sort_index
    2575       11324 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: icgf_set, isgf_set
    2576       11324 :       INTEGER, DIMENSION(:), POINTER                     :: lx, ly, lz, m, npgf
    2577       11324 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: tmp
    2578       11324 :       REAL(dp), DIMENSION(:), POINTER                    :: set_radius
    2579       11324 :       REAL(dp), DIMENSION(:, :), POINTER                 :: zet
    2580       11324 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: norm_cgf
    2581       11324 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cphi, scon, sphi
    2582             : 
    2583       11324 :       NULLIFY (set_radius, zet)
    2584             : 
    2585       11324 :       IF (sort_method == basis_sort_default) RETURN
    2586             : 
    2587             :       CALL get_gto_basis_set(gto_basis_set=basis_set, &
    2588             :                              nset=nset, &
    2589             :                              maxshell=maxshell, &
    2590             :                              maxpgf=maxpgf, &
    2591             :                              maxco=maxco, &
    2592             :                              ncgf=ncgf, &
    2593             :                              npgf=npgf, &
    2594             :                              set_radius=set_radius, &
    2595         712 :                              zet=zet)
    2596             : 
    2597        2136 :       ALLOCATE (sort_index(nset))
    2598        2136 :       ALLOCATE (tmp(nset))
    2599             :       SELECT CASE (sort_method)
    2600             :       CASE (basis_sort_zet)
    2601        4460 :          DO iset = 1, nset
    2602       13724 :             tmp(iset) = MINVAL(basis_set%zet(:npgf(iset), iset))
    2603             :          END DO
    2604             :       CASE DEFAULT
    2605         712 :          CPABORT("Request basis sort criterion not implemented.")
    2606             :       END SELECT
    2607             : 
    2608         712 :       CALL sort(tmp(1:nset), nset, sort_index)
    2609             : 
    2610         712 :       ic_max = 0
    2611         712 :       is_max = 0
    2612        4460 :       DO iset = 1, nset
    2613        3748 :          ic = 0
    2614        3748 :          is = 0
    2615       10018 :          DO ishell = 1, basis_set%nshell(iset)
    2616       21952 :             DO ico = 1, nco(basis_set%l(ishell, iset))
    2617       16394 :                ic = ic + 1
    2618        5558 :                IF (ic > ic_max) ic_max = ic
    2619             :             END DO
    2620        5558 :             lshell = basis_set%l(ishell, iset)
    2621       23944 :             DO mm = -lshell, lshell
    2622       14638 :                is = is + 1
    2623        5558 :                IF (is > is_max) is_max = is
    2624             :             END DO
    2625             :          END DO
    2626             :       END DO
    2627             : 
    2628         712 :       icgf = 0
    2629         712 :       isgf = 0
    2630        2848 :       ALLOCATE (icgf_set(nset, ic_max))
    2631       41660 :       icgf_set(:, :) = 0
    2632        2848 :       ALLOCATE (isgf_set(nset, is_max))
    2633       34030 :       isgf_set(:, :) = 0
    2634             : 
    2635        4460 :       DO iset = 1, nset
    2636        3748 :          ic = 0
    2637        3748 :          is = 0
    2638       10018 :          DO ishell = 1, basis_set%nshell(iset)
    2639       21952 :             DO ico = 1, nco(basis_set%l(ishell, iset))
    2640       16394 :                icgf = icgf + 1
    2641       16394 :                ic = ic + 1
    2642       21952 :                icgf_set(iset, ic) = icgf
    2643             :             END DO
    2644        5558 :             lshell = basis_set%l(ishell, iset)
    2645       23944 :             DO mm = -lshell, lshell
    2646       14638 :                isgf = isgf + 1
    2647       14638 :                is = is + 1
    2648       20196 :                isgf_set(iset, is) = isgf
    2649             :             END DO
    2650             :          END DO
    2651             :       END DO
    2652             : 
    2653        2136 :       ALLOCATE (cgf_symbol(SIZE(basis_set%cgf_symbol)))
    2654        2136 :       ALLOCATE (norm_cgf(SIZE(basis_set%norm_cgf)))
    2655        2136 :       ALLOCATE (lx(SIZE(basis_set%lx)))
    2656        2136 :       ALLOCATE (ly(SIZE(basis_set%ly)))
    2657        2136 :       ALLOCATE (lz(SIZE(basis_set%lz)))
    2658        2848 :       ALLOCATE (cphi(SIZE(basis_set%cphi, 1), SIZE(basis_set%cphi, 2)))
    2659      603090 :       cphi = 0.0_dp
    2660        2848 :       ALLOCATE (sphi(SIZE(basis_set%sphi, 1), SIZE(basis_set%sphi, 2)))
    2661      532390 :       sphi = 0.0_dp
    2662        2848 :       ALLOCATE (scon(SIZE(basis_set%scon, 1), SIZE(basis_set%scon, 2)))
    2663      532390 :       scon = 0.0_dp
    2664             : 
    2665        2136 :       ALLOCATE (sgf_symbol(SIZE(basis_set%sgf_symbol)))
    2666        2136 :       ALLOCATE (m(SIZE(basis_set%m)))
    2667             : 
    2668             :       icgf_new = 0
    2669             :       isgf_new = 0
    2670        4460 :       DO iset = 1, nset
    2671       37124 :          DO ic = 1, ic_max
    2672       33376 :             icgf_old = icgf_set(sort_index(iset), ic)
    2673       33376 :             IF (icgf_old == 0) CYCLE
    2674       16394 :             icgf_new = icgf_new + 1
    2675       16394 :             norm_cgf(icgf_new) = basis_set%norm_cgf(icgf_old)
    2676       16394 :             lx(icgf_new) = basis_set%lx(icgf_old)
    2677       16394 :             ly(icgf_new) = basis_set%ly(icgf_old)
    2678       16394 :             lz(icgf_new) = basis_set%lz(icgf_old)
    2679      602378 :             cphi(:, icgf_new) = basis_set%cphi(:, icgf_old)
    2680       37124 :             cgf_symbol(icgf_new) = basis_set%cgf_symbol(icgf_old)
    2681             :          END DO
    2682       31126 :          DO is = 1, is_max
    2683       26666 :             isgf_old = isgf_set(sort_index(iset), is)
    2684       26666 :             IF (isgf_old == 0) CYCLE
    2685       14638 :             isgf_new = isgf_new + 1
    2686       14638 :             m(isgf_new) = basis_set%m(isgf_old)
    2687      531678 :             sphi(:, isgf_new) = basis_set%sphi(:, isgf_old)
    2688      531678 :             scon(:, isgf_new) = basis_set%scon(:, isgf_old)
    2689       30414 :             sgf_symbol(isgf_new) = basis_set%sgf_symbol(isgf_old)
    2690             :          END DO
    2691             :       END DO
    2692             : 
    2693         712 :       DEALLOCATE (basis_set%cgf_symbol)
    2694         712 :       basis_set%cgf_symbol => cgf_symbol
    2695         712 :       DEALLOCATE (basis_set%norm_cgf)
    2696         712 :       basis_set%norm_cgf => norm_cgf
    2697         712 :       DEALLOCATE (basis_set%lx)
    2698         712 :       basis_set%lx => lx
    2699         712 :       DEALLOCATE (basis_set%ly)
    2700         712 :       basis_set%ly => ly
    2701         712 :       DEALLOCATE (basis_set%lz)
    2702         712 :       basis_set%lz => lz
    2703         712 :       DEALLOCATE (basis_set%cphi)
    2704         712 :       basis_set%cphi => cphi
    2705         712 :       DEALLOCATE (basis_set%sphi)
    2706         712 :       basis_set%sphi => sphi
    2707         712 :       DEALLOCATE (basis_set%scon)
    2708         712 :       basis_set%scon => scon
    2709             : 
    2710         712 :       DEALLOCATE (basis_set%m)
    2711         712 :       basis_set%m => m
    2712         712 :       DEALLOCATE (basis_set%sgf_symbol)
    2713         712 :       basis_set%sgf_symbol => sgf_symbol
    2714             : 
    2715        8208 :       basis_set%lmax = basis_set%lmax(sort_index)
    2716        8208 :       basis_set%lmin = basis_set%lmin(sort_index)
    2717        8208 :       basis_set%npgf = basis_set%npgf(sort_index)
    2718        8208 :       basis_set%nshell = basis_set%nshell(sort_index)
    2719        8208 :       basis_set%ncgf_set = basis_set%ncgf_set(sort_index)
    2720        8208 :       basis_set%nsgf_set = basis_set%nsgf_set(sort_index)
    2721             : 
    2722       21988 :       basis_set%n(:, :) = basis_set%n(:, sort_index)
    2723       21988 :       basis_set%l(:, :) = basis_set%l(:, sort_index)
    2724       24736 :       basis_set%zet(:, :) = basis_set%zet(:, sort_index)
    2725             : 
    2726       91816 :       basis_set%gcc(:, :, :) = basis_set%gcc(:, :, sort_index)
    2727        8208 :       basis_set%set_radius(:) = basis_set%set_radius(sort_index)
    2728       24736 :       basis_set%pgf_radius(:, :) = basis_set%pgf_radius(:, sort_index)
    2729             : 
    2730         712 :       nc = 0
    2731         712 :       ns = 0
    2732        4460 :       DO iset = 1, nset
    2733       10018 :          DO ishell = 1, basis_set%nshell(iset)
    2734        5558 :             lshell = basis_set%l(ishell, iset)
    2735        5558 :             basis_set%first_cgf(ishell, iset) = nc + 1
    2736        5558 :             nc = nc + nco(lshell)
    2737        5558 :             basis_set%last_cgf(ishell, iset) = nc
    2738        5558 :             basis_set%first_sgf(ishell, iset) = ns + 1
    2739        5558 :             ns = ns + nso(lshell)
    2740        9306 :             basis_set%last_sgf(ishell, iset) = ns
    2741             :          END DO
    2742             :       END DO
    2743             : 
    2744       12036 :    END SUBROUTINE
    2745             : 
    2746           0 : END MODULE basis_set_types

Generated by: LCOV version 1.15