LCOV - code coverage report
Current view: top level - src/aobasis - basis_set_types.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:4c33f95) Lines: 1375 1456 94.4 %
Date: 2025-01-30 06:53:08 Functions: 27 32 84.4 %

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

Generated by: LCOV version 1.15