LCOV - code coverage report
Current view: top level - src/aobasis - basis_set_types.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:7641cd9) Lines: 91.5 % 1563 1430
Test Date: 2026-05-25 07:16:39 Functions: 82.4 % 34 28

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

Generated by: LCOV version 2.0-1