LCOV - code coverage report
Current view: top level - src - qs_kind_types.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 77.1 % 1830 1411
Test Date: 2025-07-25 12:55:17 Functions: 80.8 % 26 21

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief   Define the quickstep kind type and their sub types
      10              : !> \author  Ole Schuett
      11              : !>
      12              : !> <b>Modification history:</b>
      13              : !> - 01.2002 creation [MK]
      14              : !> - 04.2002 added pao [fawzi]
      15              : !> - 09.2002 adapted for POL/KG use [GT]
      16              : !> - 02.2004 flexible normalization of basis sets [jgh]
      17              : !> - 03.2004 attach/detach routines [jgh]
      18              : !> - 10.2004 removed pao [fawzi]
      19              : !> - 08.2014 separated qs-related stuff from atomic_kind_types.F [Ole Schuett]
      20              : !> - 07.2015 new container for basis sets [jgh]
      21              : !> - 04.2021 init dft_plus_u_type [MK]
      22              : ! **************************************************************************************************
      23              : MODULE qs_kind_types
      24              :    USE atom_sgp,                        ONLY: atom_sgp_potential_type,&
      25              :                                               atom_sgp_release,&
      26              :                                               sgp_construction
      27              :    USE atom_types,                      ONLY: atom_ecppot_type,&
      28              :                                               lmat,&
      29              :                                               read_ecp_potential
      30              :    USE atom_upf,                        ONLY: atom_read_upf,&
      31              :                                               atom_release_upf,&
      32              :                                               atom_upfpot_type
      33              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      34              :                                               get_atomic_kind
      35              :    USE basis_set_container_types,       ONLY: add_basis_set_to_container,&
      36              :                                               basis_set_container_type,&
      37              :                                               get_basis_from_container,&
      38              :                                               remove_basis_from_container,&
      39              :                                               remove_basis_set_container
      40              :    USE basis_set_types,                 ONLY: &
      41              :         allocate_gto_basis_set, allocate_sto_basis_set, combine_basis_sets, &
      42              :         create_gto_from_sto_basis, deallocate_sto_basis_set, get_gto_basis_set, &
      43              :         gto_basis_set_type, init_aux_basis_set, init_orb_basis_set, read_gto_basis_set, &
      44              :         read_sto_basis_set, sto_basis_set_type, write_gto_basis_set, write_orb_basis_set
      45              :    USE cp_control_types,                ONLY: dft_control_type,&
      46              :                                               qs_control_type,&
      47              :                                               xtb_control_type
      48              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      49              :                                               cp_logger_get_default_io_unit,&
      50              :                                               cp_logger_type
      51              :    USE cp_output_handling,              ONLY: cp_p_file,&
      52              :                                               cp_print_key_finished_output,&
      53              :                                               cp_print_key_should_output,&
      54              :                                               cp_print_key_unit_nr
      55              :    USE external_potential_types,        ONLY: &
      56              :         all_potential_type, allocate_potential, deallocate_potential, get_potential, &
      57              :         gth_potential_type, init_potential, local_potential_type, read_potential, &
      58              :         set_default_all_potential, set_potential, sgp_potential_type, write_potential
      59              :    USE gapw_1c_basis_set,               ONLY: create_1c_basis
      60              :    USE input_constants,                 ONLY: &
      61              :         do_method_am1, do_method_dftb, do_method_mndo, do_method_mndod, do_method_pdg, &
      62              :         do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_pw, &
      63              :         do_method_rm1, do_method_xtb, do_qs, do_sirius, gapw_1c_large, gapw_1c_medium, &
      64              :         gapw_1c_orb, gapw_1c_small, gapw_1c_very_large
      65              :    USE input_section_types,             ONLY: section_vals_get,&
      66              :                                               section_vals_get_subs_vals,&
      67              :                                               section_vals_type,&
      68              :                                               section_vals_val_get
      69              :    USE kinds,                           ONLY: default_path_length,&
      70              :                                               default_string_length,&
      71              :                                               dp
      72              :    USE mathconstants,                   ONLY: pi
      73              :    USE message_passing,                 ONLY: mp_para_env_type
      74              :    USE orbital_pointers,                ONLY: init_orbital_pointers,&
      75              :                                               nco,&
      76              :                                               ncoset
      77              :    USE paw_proj_set_types,              ONLY: allocate_paw_proj_set,&
      78              :                                               deallocate_paw_proj_set,&
      79              :                                               get_paw_proj_set,&
      80              :                                               paw_proj_set_type,&
      81              :                                               projectors
      82              :    USE periodic_table,                  ONLY: get_ptable_info,&
      83              :                                               ptable
      84              :    USE physcon,                         ONLY: angstrom,&
      85              :                                               bohr,&
      86              :                                               evolt
      87              :    USE qs_dftb_types,                   ONLY: qs_dftb_atom_type
      88              :    USE qs_dftb_utils,                   ONLY: deallocate_dftb_atom_param,&
      89              :                                               get_dftb_atom_param,&
      90              :                                               write_dftb_atom_param
      91              :    USE qs_dispersion_types,             ONLY: qs_atom_dispersion_type
      92              :    USE qs_grid_atom,                    ONLY: allocate_grid_atom,&
      93              :                                               deallocate_grid_atom,&
      94              :                                               grid_atom_type
      95              :    USE qs_harmonics_atom,               ONLY: allocate_harmonics_atom,&
      96              :                                               deallocate_harmonics_atom,&
      97              :                                               harmonics_atom_type
      98              :    USE semi_empirical_types,            ONLY: get_se_param,&
      99              :                                               semi_empirical_create,&
     100              :                                               semi_empirical_release,&
     101              :                                               semi_empirical_type,&
     102              :                                               write_se_param
     103              :    USE semi_empirical_utils,            ONLY: init_se_param,&
     104              :                                               se_param_set_default
     105              :    USE soft_basis_set,                  ONLY: create_soft_basis
     106              :    USE string_utilities,                ONLY: uppercase
     107              :    USE xtb_parameters,                  ONLY: xtb_set_kab
     108              :    USE xtb_types,                       ONLY: deallocate_xtb_atom_param,&
     109              :                                               get_xtb_atom_param,&
     110              :                                               write_xtb_atom_param,&
     111              :                                               xtb_atom_type
     112              : #include "./base/base_uses.f90"
     113              : 
     114              :    IMPLICIT NONE
     115              : 
     116              :    PRIVATE
     117              : 
     118              :    ! Global parameters (only in this module)
     119              : 
     120              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kind_types'
     121              : 
     122              : ! **************************************************************************************************
     123              : !> \brief Input parameters for the DFT+U method
     124              : ! **************************************************************************************************
     125              :    TYPE dft_plus_u_type
     126              :       INTEGER                                :: l = -1
     127              :       INTEGER                                :: n = -1
     128              :       INTEGER                                :: max_scf = -1
     129              :       REAL(KIND=dp)                          :: eps_u_ramping = 0.0_dp
     130              :       REAL(KIND=dp)                          :: eps_scf = HUGE(0.0_dp)
     131              :       REAL(KIND=dp)                          :: u_minus_j_target = 0.0_dp
     132              :       REAL(KIND=dp)                          :: u_minus_j = 0.0_dp
     133              :       REAL(KIND=dp)                          :: u_ramping = 0.0_dp
     134              :       REAL(KIND=dp)                          :: U = 0.0_dp
     135              :       REAL(KIND=dp)                          :: J = 0.0_dp
     136              :       REAL(KIND=dp)                          :: alpha = 0.0_dp
     137              :       REAL(KIND=dp)                          :: beta = 0.0_dp
     138              :       REAL(KIND=dp)                          :: J0 = 0.0_dp
     139              :       REAL(KIND=dp)                          :: occupation = -1.0_dp
     140              :       INTEGER, DIMENSION(:), POINTER         :: orbitals => Null()
     141              :       LOGICAL                                :: init_u_ramping_each_scf = .FALSE.
     142              :       LOGICAL                                :: smear = .FALSE.
     143              :       REAL(KIND=dp), DIMENSION(:), POINTER   :: nelec => Null()
     144              :    END TYPE dft_plus_u_type
     145              : 
     146              : ! **************************************************************************************************
     147              : !> \brief Holds information about a PAO potential
     148              : ! **************************************************************************************************
     149              :    TYPE pao_potential_type
     150              :       INTEGER                                :: maxl = -1
     151              :       REAL(KIND=dp)                          :: beta = 0.0_dp
     152              :       REAL(KIND=dp)                          :: weight = 0.0_dp
     153              :       INTEGER                                :: max_projector = -1
     154              :       REAL(KIND=dp)                          :: beta_radius = HUGE(dp)
     155              :    END TYPE pao_potential_type
     156              : 
     157              : ! **************************************************************************************************
     158              : !> \brief Holds information about a PAO descriptor
     159              : ! **************************************************************************************************
     160              :    TYPE pao_descriptor_type
     161              :       REAL(KIND=dp)                          :: beta = 0.0_dp
     162              :       REAL(KIND=dp)                          :: beta_radius = HUGE(dp)
     163              :       REAL(KIND=dp)                          :: weight = 0.0_dp
     164              :       REAL(KIND=dp)                          :: screening = 0.0_dp
     165              :       REAL(KIND=dp)                          :: screening_radius = HUGE(dp)
     166              :    END TYPE pao_descriptor_type
     167              : 
     168              : ! **************************************************************************************************
     169              : !> \brief Provides all information about a quickstep kind
     170              : ! **************************************************************************************************
     171              :    TYPE qs_kind_type
     172              :       CHARACTER(LEN=default_string_length)   :: name = ""
     173              :       CHARACTER(LEN=2)                       :: element_symbol = ""
     174              :       INTEGER                                :: natom = -1
     175              :       TYPE(all_potential_type), POINTER      :: all_potential => Null()
     176              :       TYPE(local_potential_type), POINTER    :: tnadd_potential => Null()
     177              :       TYPE(gth_potential_type), POINTER      :: gth_potential => Null()
     178              :       TYPE(sgp_potential_type), POINTER      :: sgp_potential => Null()
     179              :       TYPE(semi_empirical_type), POINTER     :: se_parameter => Null()
     180              :       TYPE(qs_dftb_atom_type), POINTER       :: dftb_parameter => Null()
     181              :       TYPE(xtb_atom_type), POINTER           :: xtb_parameter => Null()
     182              :       !
     183              :       TYPE(atom_upfpot_type), POINTER        :: upf_potential => Null()
     184              :       !
     185              :       TYPE(basis_set_container_type), &
     186              :          DIMENSION(20)                       :: basis_sets = basis_set_container_type()
     187              :       ! Atomic radii
     188              :       REAL(KIND=dp)                          :: covalent_radius = 0.0_dp
     189              :       REAL(KIND=dp)                          :: vdw_radius = 0.0_dp
     190              :       ! GAPW specific data
     191              :       TYPE(paw_proj_set_type), POINTER       :: paw_proj_set => Null()
     192              :       REAL(KIND=dp)                          :: hard_radius = 0.8_dp*bohr ! for hard and soft exp
     193              :       REAL(KIND=dp)                          :: hard0_radius = 0.8_dp*bohr ! for hard exp of rho0
     194              :       REAL(KIND=dp)                          :: max_rad_local = 13.2_dp*bohr ! max GTO radius used in GAPW
     195              :       LOGICAL                                :: paw_atom = .FALSE. ! needs atomic rho1
     196              :       LOGICAL                                :: gpw_type_forced = .FALSE. ! gpw atom even if with hard exponents
     197              :       !
     198              :       LOGICAL                                :: ghost = .FALSE.
     199              :       LOGICAL                                :: floating = .FALSE.
     200              :       INTEGER                                :: lmax_dftb = -1
     201              :       REAL(KIND=dp)                          :: dudq_dftb3 = 0.0_dp
     202              :       REAL(KIND=dp)                          :: magnetization = 0.0_dp
     203              :       INTEGER, DIMENSION(:, :), POINTER      :: addel => Null()
     204              :       INTEGER, DIMENSION(:, :), POINTER      :: laddel => Null()
     205              :       INTEGER, DIMENSION(:, :), POINTER      :: naddel => Null()
     206              :       TYPE(harmonics_atom_type), POINTER     :: harmonics => Null()
     207              :       TYPE(grid_atom_type), POINTER          :: grid_atom => Null()
     208              :       INTEGER                                :: ngrid_rad = 50
     209              :       INTEGER                                :: ngrid_ang = 50
     210              :       INTEGER                                :: lmax_rho0 = 0
     211              :       INTEGER                                :: mao = -1
     212              :       INTEGER, DIMENSION(:), POINTER         :: elec_conf => Null() ! used to set up the initial atomic guess
     213              :       LOGICAL                                :: bs_occupation = .FALSE.
     214              :       TYPE(dft_plus_u_type), POINTER         :: dft_plus_u => Null()
     215              :       LOGICAL                                :: no_optimize = .TRUE.
     216              :       !
     217              :       REAL(KIND=dp), DIMENSION(:, :), POINTER :: nlcc_pot => Null()
     218              :       !
     219              :       TYPE(qs_atom_dispersion_type), POINTER :: dispersion => Null()
     220              :       REAL(KIND=dp), DIMENSION(:, :), POINTER :: reltmat => Null()
     221              :       INTEGER                                :: pao_basis_size = -1
     222              :       CHARACTER(LEN=default_path_length)     :: pao_model_file = ""
     223              :       TYPE(pao_potential_type), DIMENSION(:), POINTER :: pao_potentials => Null()
     224              :       TYPE(pao_descriptor_type), DIMENSION(:), POINTER :: pao_descriptors => Null()
     225              :    END TYPE qs_kind_type
     226              : 
     227              : ! **************************************************************************************************
     228              : !> \brief Provides a vector of pointers of type qs_kind_type
     229              : ! **************************************************************************************************
     230              :    TYPE qs_kind_p_type
     231              :       TYPE(qs_kind_type), DIMENSION(:), &
     232              :          POINTER                             :: qs_kind_set => NULL()
     233              :    END TYPE qs_kind_p_type
     234              : 
     235              :    ! Public subroutines
     236              : 
     237              :    PUBLIC :: check_qs_kind_set, &
     238              :              deallocate_qs_kind_set, &
     239              :              get_qs_kind, &
     240              :              get_qs_kind_set, &
     241              :              has_nlcc, &
     242              :              init_qs_kind_set, &
     243              :              init_gapw_basis_set, &
     244              :              init_gapw_nlcc, &
     245              :              create_qs_kind_set, &
     246              :              set_qs_kind, &
     247              :              write_qs_kind_set, &
     248              :              write_gto_basis_sets, &
     249              :              init_atom_electronic_state, set_pseudo_state
     250              : 
     251              :    ! Public data types
     252              :    PUBLIC :: qs_kind_type, pao_potential_type, pao_descriptor_type
     253              : 
     254              : CONTAINS
     255              : 
     256              : ! **************************************************************************************************
     257              : !> \brief   Destructor routine for a set of qs kinds
     258              : !> \param qs_kind_set ...
     259              : !> \date    02.01.2002
     260              : !> \author  Matthias Krack (MK)
     261              : !> \version 2.0
     262              : ! **************************************************************************************************
     263         7540 :    SUBROUTINE deallocate_qs_kind_set(qs_kind_set)
     264              : 
     265              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     266              : 
     267              :       INTEGER                                            :: ikind, nkind
     268              : 
     269         7540 :       IF (ASSOCIATED(qs_kind_set)) THEN
     270              : 
     271         7540 :          nkind = SIZE(qs_kind_set)
     272              : 
     273        22113 :          DO ikind = 1, nkind
     274        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%all_potential)) THEN
     275         6022 :                CALL deallocate_potential(qs_kind_set(ikind)%all_potential)
     276              :             END IF
     277        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%tnadd_potential)) THEN
     278           20 :                CALL deallocate_potential(qs_kind_set(ikind)%tnadd_potential)
     279              :             END IF
     280        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%gth_potential)) THEN
     281         8339 :                CALL deallocate_potential(qs_kind_set(ikind)%gth_potential)
     282              :             END IF
     283        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%sgp_potential)) THEN
     284           40 :                CALL deallocate_potential(qs_kind_set(ikind)%sgp_potential)
     285              :             END IF
     286        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%upf_potential)) THEN
     287           20 :                CALL atom_release_upf(qs_kind_set(ikind)%upf_potential)
     288           20 :                DEALLOCATE (qs_kind_set(ikind)%upf_potential)
     289              :             END IF
     290        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%se_parameter)) THEN
     291         2240 :                CALL semi_empirical_release(qs_kind_set(ikind)%se_parameter)
     292              :             END IF
     293        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%dftb_parameter)) THEN
     294          480 :                CALL deallocate_dftb_atom_param(qs_kind_set(ikind)%dftb_parameter)
     295              :             END IF
     296        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%xtb_parameter)) THEN
     297         2208 :                CALL deallocate_xtb_atom_param(qs_kind_set(ikind)%xtb_parameter)
     298              :             END IF
     299        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%paw_proj_set)) THEN
     300         1656 :                CALL deallocate_paw_proj_set(qs_kind_set(ikind)%paw_proj_set)
     301              :             END IF
     302        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%harmonics)) THEN
     303         1978 :                CALL deallocate_harmonics_atom(qs_kind_set(ikind)%harmonics)
     304              :             END IF
     305        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%grid_atom)) THEN
     306         1978 :                CALL deallocate_grid_atom(qs_kind_set(ikind)%grid_atom)
     307              :             END IF
     308        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%elec_conf)) THEN
     309        14247 :                DEALLOCATE (qs_kind_set(ikind)%elec_conf)
     310              :             END IF
     311              : 
     312        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u)) THEN
     313           32 :                IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u%orbitals)) THEN
     314            4 :                   DEALLOCATE (qs_kind_set(ikind)%dft_plus_u%orbitals)
     315              :                END IF
     316           32 :                IF (ASSOCIATED(qs_kind_set(ikind)%dft_plus_u%nelec)) THEN
     317            4 :                   DEALLOCATE (qs_kind_set(ikind)%dft_plus_u%nelec)
     318              :                END IF
     319           32 :                DEALLOCATE (qs_kind_set(ikind)%dft_plus_u)
     320              :             END IF
     321              : 
     322        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%nlcc_pot)) THEN
     323            2 :                DEALLOCATE (qs_kind_set(ikind)%nlcc_pot)
     324              :             END IF
     325              : 
     326        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%dispersion)) THEN
     327         2346 :                DEALLOCATE (qs_kind_set(ikind)%dispersion)
     328              :             END IF
     329        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%addel)) THEN
     330           62 :                DEALLOCATE (qs_kind_set(ikind)%addel)
     331              :             END IF
     332        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%naddel)) THEN
     333           62 :                DEALLOCATE (qs_kind_set(ikind)%naddel)
     334              :             END IF
     335        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%laddel)) THEN
     336           62 :                DEALLOCATE (qs_kind_set(ikind)%laddel)
     337              :             END IF
     338        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%reltmat)) THEN
     339           26 :                DEALLOCATE (qs_kind_set(ikind)%reltmat)
     340              :             END IF
     341              : 
     342        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%pao_potentials)) THEN
     343         9637 :                DEALLOCATE (qs_kind_set(ikind)%pao_potentials)
     344              :             END IF
     345        14573 :             IF (ASSOCIATED(qs_kind_set(ikind)%pao_descriptors)) THEN
     346         9637 :                DEALLOCATE (qs_kind_set(ikind)%pao_descriptors)
     347              :             END IF
     348              : 
     349        22113 :             CALL remove_basis_set_container(qs_kind_set(ikind)%basis_sets)
     350              : 
     351              :          END DO
     352         7540 :          DEALLOCATE (qs_kind_set)
     353              :       ELSE
     354              :          CALL cp_abort(__LOCATION__, &
     355              :                        "The pointer qs_kind_set is not associated and "// &
     356            0 :                        "cannot be deallocated")
     357              :       END IF
     358              : 
     359         7540 :    END SUBROUTINE deallocate_qs_kind_set
     360              : 
     361              : ! **************************************************************************************************
     362              : !> \brief Get attributes of an atomic kind.
     363              : !> \param qs_kind ...
     364              : !> \param basis_set ...
     365              : !> \param basis_type ...
     366              : !> \param ncgf ...
     367              : !> \param nsgf ...
     368              : !> \param all_potential ...
     369              : !> \param tnadd_potential ...
     370              : !> \param gth_potential ...
     371              : !> \param sgp_potential ...
     372              : !> \param upf_potential ...
     373              : !> \param se_parameter ...
     374              : !> \param dftb_parameter ...
     375              : !> \param xtb_parameter ...
     376              : !> \param dftb3_param ...
     377              : !> \param zatom ...
     378              : !> \param zeff ...
     379              : !> \param elec_conf ...
     380              : !> \param mao ...
     381              : !> \param lmax_dftb ...
     382              : !> \param alpha_core_charge ...
     383              : !> \param ccore_charge ...
     384              : !> \param core_charge ...
     385              : !> \param core_charge_radius ...
     386              : !> \param paw_proj_set ...
     387              : !> \param paw_atom ...
     388              : !> \param hard_radius ...
     389              : !> \param hard0_radius ...
     390              : !> \param max_rad_local ...
     391              : !> \param covalent_radius ...
     392              : !> \param vdw_radius ...
     393              : !> \param gpw_type_forced ...
     394              : !> \param harmonics ...
     395              : !> \param max_iso_not0 ...
     396              : !> \param max_s_harm ...
     397              : !> \param grid_atom ...
     398              : !> \param ngrid_ang ...
     399              : !> \param ngrid_rad ...
     400              : !> \param lmax_rho0 ...
     401              : !> \param dft_plus_u_atom ...
     402              : !> \param l_of_dft_plus_u ...
     403              : !> \param n_of_dft_plus_u ...
     404              : !> \param u_minus_j ...
     405              : !> \param U_of_dft_plus_u ...
     406              : !> \param J_of_dft_plus_u ...
     407              : !> \param alpha_of_dft_plus_u ...
     408              : !> \param beta_of_dft_plus_u ...
     409              : !> \param J0_of_dft_plus_u ...
     410              : !> \param occupation_of_dft_plus_u ...
     411              : !> \param dispersion ...
     412              : !> \param bs_occupation ...
     413              : !> \param magnetization ...
     414              : !> \param no_optimize ...
     415              : !> \param addel ...
     416              : !> \param laddel ...
     417              : !> \param naddel ...
     418              : !> \param orbitals ...
     419              : !> \param max_scf ...
     420              : !> \param eps_scf ...
     421              : !> \param smear ...
     422              : !> \param u_ramping ...
     423              : !> \param u_minus_j_target ...
     424              : !> \param eps_u_ramping ...
     425              : !> \param init_u_ramping_each_scf ...
     426              : !> \param reltmat ...
     427              : !> \param ghost ...
     428              : !> \param floating ...
     429              : !> \param name ...
     430              : !> \param element_symbol ...
     431              : !> \param pao_basis_size ...
     432              : !> \param pao_model_file ...
     433              : !> \param pao_potentials ...
     434              : !> \param pao_descriptors ...
     435              : !> \param nelec ...
     436              : ! **************************************************************************************************
     437     59434756 :    SUBROUTINE get_qs_kind(qs_kind, &
     438              :                           basis_set, basis_type, ncgf, nsgf, &
     439              :                           all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, &
     440              :                           se_parameter, dftb_parameter, xtb_parameter, &
     441              :                           dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, &
     442              :                           alpha_core_charge, ccore_charge, core_charge, core_charge_radius, &
     443              :                           paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, &
     444              :                           covalent_radius, vdw_radius, &
     445              :                           gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, &
     446              :                           ngrid_ang, ngrid_rad, lmax_rho0, &
     447              :                           dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, &
     448              :                           u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, &
     449              :                           alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, &
     450              :                           bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, &
     451              :                           max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, &
     452              :                           init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, &
     453              :                           pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
     454              : 
     455              :       TYPE(qs_kind_type)                                 :: qs_kind
     456              :       TYPE(gto_basis_set_type), OPTIONAL, POINTER        :: basis_set
     457              :       CHARACTER(len=*), OPTIONAL                         :: basis_type
     458              :       INTEGER, INTENT(OUT), OPTIONAL                     :: ncgf, nsgf
     459              :       TYPE(all_potential_type), OPTIONAL, POINTER        :: all_potential
     460              :       TYPE(local_potential_type), OPTIONAL, POINTER      :: tnadd_potential
     461              :       TYPE(gth_potential_type), OPTIONAL, POINTER        :: gth_potential
     462              :       TYPE(sgp_potential_type), OPTIONAL, POINTER        :: sgp_potential
     463              :       TYPE(atom_upfpot_type), OPTIONAL, POINTER          :: upf_potential
     464              :       TYPE(semi_empirical_type), OPTIONAL, POINTER       :: se_parameter
     465              :       TYPE(qs_dftb_atom_type), OPTIONAL, POINTER         :: dftb_parameter
     466              :       TYPE(xtb_atom_type), OPTIONAL, POINTER             :: xtb_parameter
     467              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: dftb3_param
     468              :       INTEGER, INTENT(OUT), OPTIONAL                     :: zatom
     469              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: zeff
     470              :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: elec_conf
     471              :       INTEGER, INTENT(OUT), OPTIONAL                     :: mao, lmax_dftb
     472              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: alpha_core_charge, ccore_charge, &
     473              :                                                             core_charge, core_charge_radius
     474              :       TYPE(paw_proj_set_type), OPTIONAL, POINTER         :: paw_proj_set
     475              :       LOGICAL, INTENT(OUT), OPTIONAL                     :: paw_atom
     476              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: hard_radius, hard0_radius, &
     477              :                                                             max_rad_local, covalent_radius, &
     478              :                                                             vdw_radius
     479              :       LOGICAL, INTENT(OUT), OPTIONAL                     :: gpw_type_forced
     480              :       TYPE(harmonics_atom_type), OPTIONAL, POINTER       :: harmonics
     481              :       INTEGER, INTENT(OUT), OPTIONAL                     :: max_iso_not0, max_s_harm
     482              :       TYPE(grid_atom_type), OPTIONAL, POINTER            :: grid_atom
     483              :       INTEGER, INTENT(OUT), OPTIONAL                     :: ngrid_ang, ngrid_rad, lmax_rho0
     484              :       LOGICAL, INTENT(OUT), OPTIONAL                     :: dft_plus_u_atom
     485              :       INTEGER, INTENT(OUT), OPTIONAL                     :: l_of_dft_plus_u, n_of_dft_plus_u
     486              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL :: u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, &
     487              :          alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u
     488              :       TYPE(qs_atom_dispersion_type), OPTIONAL, POINTER   :: dispersion
     489              :       LOGICAL, INTENT(OUT), OPTIONAL                     :: bs_occupation
     490              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: magnetization
     491              :       LOGICAL, INTENT(OUT), OPTIONAL                     :: no_optimize
     492              :       INTEGER, DIMENSION(:, :), OPTIONAL, POINTER        :: addel, laddel, naddel
     493              :       INTEGER, DIMENSION(:), OPTIONAL, POINTER           :: orbitals
     494              :       INTEGER, OPTIONAL                                  :: max_scf
     495              :       REAL(KIND=dp), OPTIONAL                            :: eps_scf
     496              :       LOGICAL, OPTIONAL                                  :: smear
     497              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: u_ramping, u_minus_j_target, &
     498              :                                                             eps_u_ramping
     499              :       LOGICAL, OPTIONAL                                  :: init_u_ramping_each_scf
     500              :       REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: reltmat
     501              :       LOGICAL, OPTIONAL                                  :: ghost, floating
     502              :       CHARACTER(LEN=default_string_length), &
     503              :          INTENT(OUT), OPTIONAL                           :: name
     504              :       CHARACTER(LEN=2), INTENT(OUT), OPTIONAL            :: element_symbol
     505              :       INTEGER, INTENT(OUT), OPTIONAL                     :: pao_basis_size
     506              :       CHARACTER(LEN=default_path_length), INTENT(OUT), &
     507              :          OPTIONAL                                        :: pao_model_file
     508              :       TYPE(pao_potential_type), DIMENSION(:), OPTIONAL, &
     509              :          POINTER                                         :: pao_potentials
     510              :       TYPE(pao_descriptor_type), DIMENSION(:), &
     511              :          OPTIONAL, POINTER                               :: pao_descriptors
     512              :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: nelec
     513              : 
     514              :       CHARACTER(LEN=default_string_length)               :: my_basis_type
     515              :       INTEGER                                            :: l
     516              :       LOGICAL                                            :: found
     517              :       TYPE(gto_basis_set_type), POINTER                  :: tmp_basis_set
     518              : 
     519              :       ! Retrieve basis set from the kind container
     520     59434756 :       IF (PRESENT(basis_type)) THEN
     521      9264141 :          my_basis_type = basis_type
     522              :       ELSE
     523     50170615 :          my_basis_type = "ORB"
     524              :       END IF
     525              : 
     526     59434756 :       IF (PRESENT(basis_set)) THEN
     527              :          CALL get_basis_from_container(qs_kind%basis_sets, basis_set=basis_set, &
     528      9901770 :                                        basis_type=my_basis_type)
     529              :       END IF
     530              : 
     531     59434756 :       IF (PRESENT(ncgf)) THEN
     532              :          CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
     533          960 :                                        basis_type=my_basis_type)
     534          960 :          IF (ASSOCIATED(tmp_basis_set)) THEN
     535          960 :             CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, ncgf=ncgf)
     536            0 :          ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
     537            0 :             l = qs_kind%dftb_parameter%lmax
     538            0 :             ncgf = ((l + 1)*(l + 2)*(l + 3))/6
     539              :          ELSE
     540            0 :             ncgf = 0
     541              :          END IF
     542              :       END IF
     543              : 
     544     59434756 :       IF (PRESENT(nsgf)) THEN
     545              :          CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
     546       266941 :                                        basis_type=my_basis_type)
     547       266941 :          IF (ASSOCIATED(tmp_basis_set)) THEN
     548       165127 :             CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nsgf=nsgf)
     549       101814 :          ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
     550       101810 :             nsgf = qs_kind%dftb_parameter%natorb
     551              :          ELSE
     552            4 :             nsgf = 0
     553              :          END IF
     554              :       END IF
     555              : 
     556     59434756 :       IF (PRESENT(all_potential)) all_potential => qs_kind%all_potential
     557     59434756 :       IF (PRESENT(tnadd_potential)) tnadd_potential => qs_kind%tnadd_potential
     558     59434756 :       IF (PRESENT(gth_potential)) gth_potential => qs_kind%gth_potential
     559     59434756 :       IF (PRESENT(sgp_potential)) sgp_potential => qs_kind%sgp_potential
     560     59434756 :       IF (PRESENT(upf_potential)) upf_potential => qs_kind%upf_potential
     561     59434756 :       IF (PRESENT(se_parameter)) se_parameter => qs_kind%se_parameter
     562     59434756 :       IF (PRESENT(dftb_parameter)) dftb_parameter => qs_kind%dftb_parameter
     563     59434756 :       IF (PRESENT(xtb_parameter)) xtb_parameter => qs_kind%xtb_parameter
     564     59434756 :       IF (PRESENT(element_symbol)) element_symbol = qs_kind%element_symbol
     565     59434756 :       IF (PRESENT(name)) name = qs_kind%name
     566     59434756 :       IF (PRESENT(dftb3_param)) dftb3_param = qs_kind%dudq_dftb3
     567     59434756 :       IF (PRESENT(elec_conf)) elec_conf => qs_kind%elec_conf
     568     59434756 :       IF (PRESENT(alpha_core_charge)) THEN
     569       207081 :          IF (ASSOCIATED(qs_kind%all_potential)) THEN
     570              :             CALL get_potential(potential=qs_kind%all_potential, &
     571        52298 :                                alpha_core_charge=alpha_core_charge)
     572       154783 :          ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
     573              :             CALL get_potential(potential=qs_kind%gth_potential, &
     574       152911 :                                alpha_core_charge=alpha_core_charge)
     575         1872 :          ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
     576              :             CALL get_potential(potential=qs_kind%sgp_potential, &
     577          456 :                                alpha_core_charge=alpha_core_charge)
     578              :          ELSE
     579         1416 :             alpha_core_charge = 1.0_dp
     580              :          END IF
     581              :       END IF
     582     59434756 :       IF (PRESENT(ccore_charge)) THEN
     583        84781 :          IF (ASSOCIATED(qs_kind%all_potential)) THEN
     584              :             CALL get_potential(potential=qs_kind%all_potential, &
     585        12204 :                                ccore_charge=ccore_charge)
     586        72577 :          ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
     587              :             CALL get_potential(potential=qs_kind%gth_potential, &
     588        71529 :                                ccore_charge=ccore_charge)
     589         1048 :          ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
     590              :             CALL get_potential(potential=qs_kind%sgp_potential, &
     591          234 :                                ccore_charge=ccore_charge)
     592          814 :          ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
     593            0 :             CPABORT("UPF CCORE CHARGE RADIUS NOT AVAILABLE")
     594              :          ELSE
     595          814 :             ccore_charge = 0.0_dp
     596              :          END IF
     597              :       END IF
     598     59434756 :       IF (PRESENT(core_charge_radius)) THEN
     599        85741 :          IF (ASSOCIATED(qs_kind%all_potential)) THEN
     600              :             CALL get_potential(potential=qs_kind%all_potential, &
     601        36804 :                                core_charge_radius=core_charge_radius)
     602        48937 :          ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
     603              :             CALL get_potential(potential=qs_kind%gth_potential, &
     604        48417 :                                core_charge_radius=core_charge_radius)
     605          520 :          ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
     606              :             CALL get_potential(potential=qs_kind%sgp_potential, &
     607          130 :                                core_charge_radius=core_charge_radius)
     608          390 :          ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
     609            0 :             CPABORT("UPF CORE CHARGE RADIUS NOT AVAILABLE")
     610              :          ELSE
     611          390 :             core_charge_radius = 0.0_dp
     612              :          END IF
     613              :       END IF
     614     59434756 :       IF (PRESENT(core_charge)) THEN
     615        39741 :          IF (ASSOCIATED(qs_kind%all_potential)) THEN
     616              :             CALL get_potential(potential=qs_kind%all_potential, &
     617          367 :                                zeff=core_charge)
     618        39374 :          ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
     619              :             CALL get_potential(potential=qs_kind%gth_potential, &
     620        39374 :                                zeff=core_charge)
     621            0 :          ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
     622              :             CALL get_potential(potential=qs_kind%sgp_potential, &
     623            0 :                                zeff=core_charge)
     624            0 :          ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
     625            0 :             CPABORT("UPF CORE CHARGE NOT AVAILABLE")
     626              :          ELSE
     627            0 :             core_charge = 0.0_dp
     628              :          END IF
     629              :       END IF
     630              : 
     631     59434756 :       IF (PRESENT(zatom)) THEN
     632              :          ! Retrieve information on element
     633       225238 :          CALL get_ptable_info(qs_kind%element_symbol, ielement=zatom, found=found)
     634       225238 :          CPASSERT(found)
     635              :       END IF
     636              : 
     637     59434756 :       IF (PRESENT(zeff)) THEN
     638       222992 :          IF (ASSOCIATED(qs_kind%all_potential)) THEN
     639        56571 :             CALL get_potential(potential=qs_kind%all_potential, zeff=zeff)
     640       166421 :          ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
     641       164975 :             CALL get_potential(potential=qs_kind%gth_potential, zeff=zeff)
     642         1446 :          ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
     643          400 :             CALL get_potential(potential=qs_kind%sgp_potential, zeff=zeff)
     644         1046 :          ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
     645           65 :             zeff = qs_kind%upf_potential%zion
     646              :          ELSE
     647          981 :             zeff = 0.0_dp
     648              :          END IF
     649              :       END IF
     650              : 
     651     59434756 :       IF (PRESENT(covalent_radius)) covalent_radius = qs_kind%covalent_radius
     652     59434756 :       IF (PRESENT(vdw_radius)) vdw_radius = qs_kind%vdw_radius
     653              : 
     654     59434756 :       IF (PRESENT(paw_proj_set)) paw_proj_set => qs_kind%paw_proj_set
     655     59434756 :       IF (PRESENT(paw_atom)) paw_atom = qs_kind%paw_atom
     656     59434756 :       IF (PRESENT(gpw_type_forced)) gpw_type_forced = qs_kind%gpw_type_forced
     657     59434756 :       IF (PRESENT(hard_radius)) hard_radius = qs_kind%hard_radius
     658     59434756 :       IF (PRESENT(hard0_radius)) hard0_radius = qs_kind%hard0_radius
     659     59434756 :       IF (PRESENT(max_rad_local)) max_rad_local = qs_kind%max_rad_local
     660     59434756 :       IF (PRESENT(harmonics)) harmonics => qs_kind%harmonics
     661     59434756 :       IF (PRESENT(max_s_harm)) THEN
     662      8443742 :          IF (ASSOCIATED(qs_kind%harmonics)) THEN
     663       286308 :             max_s_harm = qs_kind%harmonics%max_s_harm
     664              :          ELSE
     665      8157434 :             max_s_harm = 0
     666              :          END IF
     667              :       END IF
     668     59434756 :       IF (PRESENT(max_iso_not0)) THEN
     669      8475108 :          IF (ASSOCIATED(qs_kind%harmonics)) THEN
     670       317674 :             max_iso_not0 = qs_kind%harmonics%max_iso_not0
     671              :          ELSE
     672      8157434 :             max_iso_not0 = 0
     673              :          END IF
     674              :       END IF
     675     59434756 :       IF (PRESENT(grid_atom)) grid_atom => qs_kind%grid_atom
     676     59434756 :       IF (PRESENT(ngrid_ang)) ngrid_ang = qs_kind%ngrid_ang
     677     59434756 :       IF (PRESENT(ngrid_rad)) ngrid_rad = qs_kind%ngrid_rad
     678     59434756 :       IF (PRESENT(lmax_rho0)) lmax_rho0 = qs_kind%lmax_rho0
     679     59434756 :       IF (PRESENT(ghost)) ghost = qs_kind%ghost
     680     59434756 :       IF (PRESENT(floating)) floating = qs_kind%floating
     681     59434756 :       IF (PRESENT(dft_plus_u_atom)) dft_plus_u_atom = ASSOCIATED(qs_kind%dft_plus_u)
     682     59434756 :       IF (PRESENT(l_of_dft_plus_u)) THEN
     683         5106 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     684         2542 :             l_of_dft_plus_u = qs_kind%dft_plus_u%l
     685              :          ELSE
     686         2564 :             l_of_dft_plus_u = -1
     687              :          END IF
     688              :       END IF
     689     59434756 :       IF (PRESENT(n_of_dft_plus_u)) THEN
     690           22 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     691            0 :             n_of_dft_plus_u = qs_kind%dft_plus_u%n
     692              :          ELSE
     693           22 :             n_of_dft_plus_u = -1
     694              :          END IF
     695              :       END IF
     696     59434756 :       IF (PRESENT(u_minus_j)) THEN
     697         5084 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     698         2542 :             u_minus_j = qs_kind%dft_plus_u%u_minus_j
     699              :          ELSE
     700         2542 :             u_minus_j = 0.0_dp
     701              :          END IF
     702              :       END IF
     703     59434756 :       IF (PRESENT(u_minus_j_target)) THEN
     704         5106 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     705         2542 :             u_minus_j_target = qs_kind%dft_plus_u%u_minus_j_target
     706              :          ELSE
     707         2564 :             u_minus_j_target = 0.0_dp
     708              :          END IF
     709              :       END IF
     710     59434756 :       IF (PRESENT(U_of_dft_plus_u)) THEN
     711           22 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     712            0 :             U_of_dft_plus_u = qs_kind%dft_plus_u%U
     713              :          ELSE
     714           22 :             U_of_dft_plus_u = 0.0_dp
     715              :          END IF
     716              :       END IF
     717     59434756 :       IF (PRESENT(J_of_dft_plus_u)) THEN
     718           22 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     719            0 :             J_of_dft_plus_u = qs_kind%dft_plus_u%J
     720              :          ELSE
     721           22 :             J_of_dft_plus_u = 0.0_dp
     722              :          END IF
     723              :       END IF
     724     59434756 :       IF (PRESENT(alpha_of_dft_plus_u)) THEN
     725           22 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     726            0 :             alpha_of_dft_plus_u = qs_kind%dft_plus_u%alpha
     727              :          ELSE
     728           22 :             alpha_of_dft_plus_u = 0.0_dp
     729              :          END IF
     730              :       END IF
     731     59434756 :       IF (PRESENT(beta_of_dft_plus_u)) THEN
     732           22 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     733            0 :             beta_of_dft_plus_u = qs_kind%dft_plus_u%beta
     734              :          ELSE
     735           22 :             beta_of_dft_plus_u = 0.0_dp
     736              :          END IF
     737              :       END IF
     738     59434756 :       IF (PRESENT(J0_of_dft_plus_u)) THEN
     739           22 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     740            0 :             J0_of_dft_plus_u = qs_kind%dft_plus_u%J0
     741              :          ELSE
     742           22 :             J0_of_dft_plus_u = 0.0_dp
     743              :          END IF
     744              :       END IF
     745     59434756 :       IF (PRESENT(occupation_of_dft_plus_u)) THEN
     746           22 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     747            0 :             occupation_of_dft_plus_u = qs_kind%dft_plus_u%occupation
     748              :          ELSE
     749           22 :             occupation_of_dft_plus_u = -1.0_dp
     750              :          END IF
     751              :       END IF
     752              : 
     753     59434756 :       IF (PRESENT(init_u_ramping_each_scf)) THEN
     754          160 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     755           80 :             init_u_ramping_each_scf = qs_kind%dft_plus_u%init_u_ramping_each_scf
     756              :          ELSE
     757           80 :             init_u_ramping_each_scf = .FALSE.
     758              :          END IF
     759              :       END IF
     760     59434756 :       IF (PRESENT(u_ramping)) THEN
     761         5244 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     762         2622 :             u_ramping = qs_kind%dft_plus_u%u_ramping
     763              :          ELSE
     764         2622 :             u_ramping = 0.0_dp
     765              :          END IF
     766              :       END IF
     767     59434756 :       IF (PRESENT(eps_u_ramping)) THEN
     768         5084 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     769         2542 :             eps_u_ramping = qs_kind%dft_plus_u%eps_u_ramping
     770              :          ELSE
     771         2542 :             eps_u_ramping = 1.0E-5_dp
     772              :          END IF
     773              :       END IF
     774     59434756 :       IF (PRESENT(nelec)) THEN
     775         3840 :          NULLIFY (nelec)
     776         3840 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     777         1920 :             IF (ASSOCIATED(qs_kind%dft_plus_u%nelec)) THEN
     778            0 :                nelec => qs_kind%dft_plus_u%nelec
     779              :             END IF
     780              :          END IF
     781              :       END IF
     782     59434756 :       IF (PRESENT(orbitals)) THEN
     783         4160 :          NULLIFY (orbitals)
     784         4160 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     785         2080 :             IF (ASSOCIATED(qs_kind%dft_plus_u%orbitals)) THEN
     786          136 :                orbitals => qs_kind%dft_plus_u%orbitals
     787              :             END IF
     788              :          END IF
     789              :       END IF
     790     59434756 :       IF (PRESENT(eps_scf)) THEN
     791         4160 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     792         2080 :             eps_scf = qs_kind%dft_plus_u%eps_scf
     793              :          ELSE
     794         2080 :             eps_scf = 1.0E30_dp
     795              :          END IF
     796              :       END IF
     797     59434756 :       IF (PRESENT(max_scf)) THEN
     798         4160 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     799         2080 :             max_scf = qs_kind%dft_plus_u%max_scf
     800              :          ELSE
     801         2080 :             max_scf = -1
     802              :          END IF
     803              :       END IF
     804     59434756 :       IF (PRESENT(smear)) THEN
     805         4160 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
     806         2080 :             smear = qs_kind%dft_plus_u%smear
     807              :          ELSE
     808         2080 :             smear = .FALSE.
     809              :          END IF
     810              :       END IF
     811     59434756 :       IF (PRESENT(dispersion)) dispersion => qs_kind%dispersion
     812     59434756 :       IF (PRESENT(bs_occupation)) bs_occupation = qs_kind%bs_occupation
     813     59434756 :       IF (PRESENT(addel)) addel => qs_kind%addel
     814     59434756 :       IF (PRESENT(laddel)) laddel => qs_kind%laddel
     815     59434756 :       IF (PRESENT(naddel)) naddel => qs_kind%naddel
     816              : 
     817     59434756 :       IF (PRESENT(magnetization)) magnetization = qs_kind%magnetization
     818              : 
     819     59434756 :       IF (PRESENT(no_optimize)) no_optimize = qs_kind%no_optimize
     820              : 
     821     59434756 :       IF (PRESENT(reltmat)) reltmat => qs_kind%reltmat
     822              : 
     823     59434756 :       IF (PRESENT(mao)) mao = qs_kind%mao
     824              : 
     825     59434756 :       IF (PRESENT(lmax_dftb)) lmax_dftb = qs_kind%lmax_dftb
     826              : 
     827     59434756 :       IF (PRESENT(pao_basis_size)) pao_basis_size = qs_kind%pao_basis_size
     828     59434756 :       IF (PRESENT(pao_model_file)) pao_model_file = qs_kind%pao_model_file
     829     59434756 :       IF (PRESENT(pao_potentials)) pao_potentials => qs_kind%pao_potentials
     830     59434756 :       IF (PRESENT(pao_descriptors)) pao_descriptors => qs_kind%pao_descriptors
     831     59434756 :    END SUBROUTINE get_qs_kind
     832              : 
     833              : ! **************************************************************************************************
     834              : !> \brief Get attributes of an atomic kind set.
     835              : !> \param qs_kind_set ...
     836              : !> \param all_potential_present ...
     837              : !> \param tnadd_potential_present ...
     838              : !> \param gth_potential_present ...
     839              : !> \param sgp_potential_present ...
     840              : !> \param paw_atom_present ...
     841              : !> \param dft_plus_u_atom_present ...
     842              : !> \param maxcgf ...
     843              : !> \param maxsgf ...
     844              : !> \param maxco ...
     845              : !> \param maxco_proj ...
     846              : !> \param maxgtops ...
     847              : !> \param maxlgto ...
     848              : !> \param maxlprj ...
     849              : !> \param maxnset ...
     850              : !> \param maxsgf_set ...
     851              : !> \param ncgf ...
     852              : !> \param npgf ...
     853              : !> \param nset ...
     854              : !> \param nsgf ...
     855              : !> \param nshell ...
     856              : !> \param maxpol ...
     857              : !> \param maxlppl ...
     858              : !> \param maxlppnl ...
     859              : !> \param maxppnl ...
     860              : !> \param nelectron ...
     861              : !> \param maxder ...
     862              : !> \param max_ngrid_rad ...
     863              : !> \param max_sph_harm ...
     864              : !> \param maxg_iso_not0 ...
     865              : !> \param lmax_rho0 ...
     866              : !> \param basis_rcut ...
     867              : !> \param basis_type ...
     868              : !> \param total_zeff_corr ... [SGh]
     869              : !> \param npgf_seg total number of primitive GTOs in "segmented contraction format"
     870              : ! **************************************************************************************************
     871      4018948 :    SUBROUTINE get_qs_kind_set(qs_kind_set, &
     872              :                               all_potential_present, tnadd_potential_present, gth_potential_present, &
     873              :                               sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, &
     874              :                               maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, &
     875              :                               ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, &
     876              :                               nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, &
     877              :                               basis_rcut, &
     878              :                               basis_type, total_zeff_corr, npgf_seg)
     879              : 
     880              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     881              :       LOGICAL, INTENT(OUT), OPTIONAL :: all_potential_present, tnadd_potential_present, &
     882              :          gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present
     883              :       INTEGER, INTENT(OUT), OPTIONAL :: maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, &
     884              :          maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, &
     885              :          maxppnl, nelectron
     886              :       INTEGER, INTENT(IN), OPTIONAL                      :: maxder
     887              :       INTEGER, INTENT(OUT), OPTIONAL                     :: max_ngrid_rad, max_sph_harm, &
     888              :                                                             maxg_iso_not0, lmax_rho0
     889              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: basis_rcut
     890              :       CHARACTER(len=*), OPTIONAL                         :: basis_type
     891              :       REAL(KIND=dp), INTENT(OUT), OPTIONAL               :: total_zeff_corr
     892              :       INTEGER, INTENT(OUT), OPTIONAL                     :: npgf_seg
     893              : 
     894              :       CHARACTER(len=default_string_length)               :: my_basis_type
     895              :       INTEGER                                            :: ikind, imax, lmax_rho0_kind, &
     896              :                                                             max_iso_not0, max_s_harm, n, &
     897              :                                                             ngrid_rad, nkind, nrloc(10), &
     898              :                                                             nrpot(1:15, 0:10)
     899              :       LOGICAL                                            :: dft_plus_u_atom, ecp_semi_local, paw_atom
     900              :       REAL(KIND=dp)                                      :: brcut, zeff, zeff_correction
     901              :       TYPE(all_potential_type), POINTER                  :: all_potential
     902              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
     903              :       TYPE(gto_basis_set_type), POINTER                  :: tmp_basis_set
     904              :       TYPE(local_potential_type), POINTER                :: tnadd_potential
     905              :       TYPE(paw_proj_set_type), POINTER                   :: paw_proj_set
     906              :       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
     907              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
     908              :       TYPE(sgp_potential_type), POINTER                  :: sgp_potential
     909              : 
     910      4018948 :       IF (PRESENT(basis_type)) THEN
     911      3729771 :          my_basis_type = basis_type
     912              :       ELSE
     913       289177 :          my_basis_type = "ORB"
     914              :       END IF
     915              : 
     916      4018948 :       IF (ASSOCIATED(qs_kind_set)) THEN
     917              : 
     918      4018948 :          IF (PRESENT(maxcgf)) maxcgf = 0
     919      4018948 :          IF (PRESENT(maxco)) maxco = 0
     920      4018948 :          IF (PRESENT(maxco_proj)) maxco_proj = 0
     921      4018948 :          IF (PRESENT(maxg_iso_not0)) maxg_iso_not0 = 0
     922      4018948 :          IF (PRESENT(maxgtops)) maxgtops = 0
     923      4018948 :          IF (PRESENT(maxlgto)) maxlgto = -1
     924      4018948 :          IF (PRESENT(maxlppl)) maxlppl = -1
     925      4018948 :          IF (PRESENT(maxlppnl)) maxlppnl = -1
     926      4018948 :          IF (PRESENT(maxpol)) maxpol = -1
     927      4018948 :          IF (PRESENT(maxlprj)) maxlprj = -1
     928      4018948 :          IF (PRESENT(maxnset)) maxnset = 0
     929      4018948 :          IF (PRESENT(maxppnl)) maxppnl = 0
     930      4018948 :          IF (PRESENT(maxsgf)) maxsgf = 0
     931      4018948 :          IF (PRESENT(maxsgf_set)) maxsgf_set = 0
     932      4018948 :          IF (PRESENT(ncgf)) ncgf = 0
     933      4018948 :          IF (PRESENT(nelectron)) nelectron = 0
     934      4018948 :          IF (PRESENT(npgf)) npgf = 0
     935      4018948 :          IF (PRESENT(nset)) nset = 0
     936      4018948 :          IF (PRESENT(nsgf)) nsgf = 0
     937      4018948 :          IF (PRESENT(nshell)) nshell = 0
     938      4018948 :          IF (PRESENT(all_potential_present)) all_potential_present = .FALSE.
     939      4018948 :          IF (PRESENT(tnadd_potential_present)) tnadd_potential_present = .FALSE.
     940      4018948 :          IF (PRESENT(gth_potential_present)) gth_potential_present = .FALSE.
     941      4018948 :          IF (PRESENT(sgp_potential_present)) sgp_potential_present = .FALSE.
     942      4018948 :          IF (PRESENT(paw_atom_present)) paw_atom_present = .FALSE.
     943      4018948 :          IF (PRESENT(max_ngrid_rad)) max_ngrid_rad = 0
     944      4018948 :          IF (PRESENT(max_sph_harm)) max_sph_harm = 0
     945      4018948 :          IF (PRESENT(lmax_rho0)) lmax_rho0 = 0
     946      4018948 :          IF (PRESENT(basis_rcut)) basis_rcut = 0.0_dp
     947      4018948 :          IF (PRESENT(total_zeff_corr)) total_zeff_corr = 0.0_dp
     948      4018948 :          IF (PRESENT(npgf_seg)) npgf_seg = 0
     949              : 
     950      4018948 :          nkind = SIZE(qs_kind_set)
     951     12462690 :          DO ikind = 1, nkind
     952      8443742 :             qs_kind => qs_kind_set(ikind)
     953              :             CALL get_qs_kind(qs_kind=qs_kind, &
     954              :                              all_potential=all_potential, &
     955              :                              tnadd_potential=tnadd_potential, &
     956              :                              gth_potential=gth_potential, &
     957              :                              sgp_potential=sgp_potential, &
     958              :                              paw_proj_set=paw_proj_set, &
     959              :                              dftb_parameter=dftb_parameter, &
     960              :                              ngrid_rad=ngrid_rad, &
     961              :                              max_s_harm=max_s_harm, &
     962              :                              max_iso_not0=max_iso_not0, &
     963              :                              paw_atom=paw_atom, &
     964              :                              dft_plus_u_atom=dft_plus_u_atom, &
     965      8443742 :                              lmax_rho0=lmax_rho0_kind)
     966              : 
     967      8443742 :             IF (PRESENT(maxlppl) .AND. ASSOCIATED(gth_potential)) THEN
     968        43527 :                CALL get_potential(potential=gth_potential, nexp_ppl=n)
     969        43527 :                maxlppl = MAX(maxlppl, 2*(n - 1))
     970         9593 :             ELSEIF (PRESENT(maxlppl) .AND. ASSOCIATED(sgp_potential)) THEN
     971          132 :                CALL get_potential(potential=sgp_potential, nrloc=nrloc, ecp_semi_local=ecp_semi_local)
     972         1452 :                n = MAXVAL(nrloc) - 2
     973          132 :                maxlppl = MAX(maxlppl, 2*(n - 1))
     974          132 :                IF (ecp_semi_local) THEN
     975          102 :                   CALL get_potential(potential=sgp_potential, sl_lmax=imax, nrpot=nrpot)
     976        18054 :                   n = MAXVAL(nrpot) - 2
     977          102 :                   n = 2*(n - 1) + imax
     978          102 :                   maxlppl = MAX(maxlppl, n)
     979              :                END IF
     980              :             END IF
     981              : 
     982      8443742 :             IF (PRESENT(maxlppnl) .AND. ASSOCIATED(gth_potential)) THEN
     983        40554 :                CALL get_potential(potential=gth_potential, lprj_ppnl_max=imax)
     984        40554 :                maxlppnl = MAX(maxlppnl, imax)
     985         9509 :             ELSEIF (PRESENT(maxlppnl) .AND. ASSOCIATED(sgp_potential)) THEN
     986           70 :                CALL get_potential(potential=sgp_potential, lmax=imax)
     987           70 :                maxlppnl = MAX(maxlppnl, imax)
     988              :             END IF
     989              : 
     990      8443742 :             IF (PRESENT(maxpol) .AND. ASSOCIATED(tnadd_potential)) THEN
     991           66 :                CALL get_potential(potential=tnadd_potential, npol=n)
     992           66 :                maxpol = MAX(maxpol, 2*(n - 1))
     993              :             END IF
     994              : 
     995      8443742 :             IF (PRESENT(maxco_proj) .AND. ASSOCIATED(paw_proj_set)) THEN
     996         4278 :                CALL get_paw_proj_set(paw_proj_set=paw_proj_set, ncgauprj=imax)
     997         4278 :                maxco_proj = MAX(maxco_proj, imax)
     998              :             END IF
     999              : 
    1000      8443742 :             IF (PRESENT(maxlprj) .AND. ASSOCIATED(paw_proj_set)) THEN
    1001         4278 :                CALL get_paw_proj_set(paw_proj_set=paw_proj_set, maxl=imax)
    1002         4278 :                maxlprj = MAX(maxlprj, imax)
    1003              :             END IF
    1004              : 
    1005      8443742 :             IF (PRESENT(maxppnl) .AND. ASSOCIATED(gth_potential)) THEN
    1006        28180 :                CALL get_potential(potential=gth_potential, nppnl=imax)
    1007        28180 :                maxppnl = MAX(maxppnl, imax)
    1008          240 :             ELSEIF (PRESENT(maxppnl) .AND. ASSOCIATED(sgp_potential)) THEN
    1009           10 :                CALL get_potential(potential=sgp_potential, nppnl=imax)
    1010           10 :                maxppnl = MAX(maxppnl, imax)
    1011              :             END IF
    1012              : 
    1013              :             CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
    1014      8443742 :                                           basis_type=my_basis_type)
    1015              : 
    1016      8443742 :             IF (PRESENT(maxcgf)) THEN
    1017            0 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1018            0 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, ncgf=imax)
    1019            0 :                   maxcgf = MAX(maxcgf, imax)
    1020            0 :                ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
    1021            0 :                   CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=imax)
    1022            0 :                   imax = ((imax + 1)*(imax + 2)*(imax + 3))/6
    1023            0 :                   maxcgf = MAX(maxcgf, imax)
    1024              :                END IF
    1025              :             END IF
    1026              : 
    1027      8443742 :             IF (PRESENT(maxco)) THEN
    1028      7711441 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1029      7711433 :                   IF (PRESENT(maxder)) THEN
    1030              :                      CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, &
    1031            0 :                                             maxco=imax, maxder=maxder)
    1032              :                   ELSE
    1033      7711433 :                      CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxco=imax)
    1034              :                   END IF
    1035      7711433 :                   maxco = MAX(maxco, imax)
    1036              :                END IF
    1037      7711441 :                IF (ASSOCIATED(gth_potential)) THEN
    1038       652823 :                   CALL get_potential(potential=gth_potential, lprj_ppnl_max=imax)
    1039       652823 :                   maxco = MAX(maxco, ncoset(imax))
    1040              :                END IF
    1041      7711441 :                IF (ASSOCIATED(sgp_potential)) THEN
    1042          992 :                   CALL get_potential(potential=sgp_potential, lmax=imax)
    1043          992 :                   maxco = MAX(maxco, ncoset(imax))
    1044          992 :                   CALL get_potential(potential=sgp_potential, sl_lmax=imax)
    1045          992 :                   maxco = MAX(maxco, ncoset(imax))
    1046              :                END IF
    1047              :             END IF
    1048              : 
    1049      8443742 :             IF (PRESENT(maxgtops)) THEN
    1050        98888 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1051        98888 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxso=imax, nset=n)
    1052        98888 :                   maxgtops = MAX(maxgtops, n*imax)
    1053              :                END IF
    1054              :             END IF
    1055              : 
    1056      8443742 :             IF (PRESENT(maxlgto)) THEN
    1057      7301624 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1058      7278669 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxl=imax)
    1059      7278669 :                   maxlgto = MAX(maxlgto, imax)
    1060        22955 :                ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
    1061         2516 :                   CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=imax)
    1062         2516 :                   maxlgto = MAX(maxlgto, imax)
    1063              :                END IF
    1064              :             END IF
    1065              : 
    1066      8443742 :             IF (PRESENT(maxnset)) THEN
    1067        75465 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1068        75465 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nset=n)
    1069        75465 :                   maxnset = MAX(maxnset, n)
    1070              :                END IF
    1071              :             END IF
    1072              : 
    1073      8443742 :             IF (PRESENT(maxsgf)) THEN
    1074      7467570 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1075      7467546 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nsgf=imax)
    1076      7467546 :                   maxsgf = MAX(maxsgf, imax)
    1077              :                END IF
    1078              :             END IF
    1079              : 
    1080      8443742 :             IF (PRESENT(maxsgf_set)) THEN
    1081       446103 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1082       446103 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, maxsgf_set=imax)
    1083       446103 :                   maxsgf_set = MAX(maxsgf_set, imax)
    1084              :                END IF
    1085              :             END IF
    1086              : 
    1087      8443742 :             IF (PRESENT(ncgf)) THEN
    1088        36443 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1089        15031 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, ncgf=n)
    1090        15031 :                   ncgf = ncgf + n*qs_kind_set(ikind)%natom
    1091        21412 :                ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
    1092          987 :                   CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=imax)
    1093          987 :                   n = ((imax + 1)*(imax + 2)*(imax + 3))/6
    1094          987 :                   ncgf = ncgf + n*qs_kind_set(ikind)%natom
    1095              :                END IF
    1096              :             END IF
    1097              : 
    1098      8443742 :             IF (PRESENT(npgf)) THEN
    1099        28800 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1100         7415 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, npgf_sum=n)
    1101         7415 :                   npgf = npgf + n*qs_kind_set(ikind)%natom
    1102              :                END IF
    1103              :             END IF
    1104              : 
    1105      8443742 :             IF (PRESENT(nset)) THEN
    1106        28800 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1107         7415 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nset=n)
    1108         7415 :                   nset = nset + n*qs_kind_set(ikind)%natom
    1109              :                END IF
    1110              :             END IF
    1111              : 
    1112      8443742 :             IF (PRESENT(nsgf)) THEN
    1113       104756 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1114        69474 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nsgf=n)
    1115        69474 :                   nsgf = nsgf + n*qs_kind_set(ikind)%natom
    1116        35282 :                ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
    1117        14855 :                   CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, natorb=n)
    1118        14855 :                   nsgf = nsgf + n*qs_kind_set(ikind)%natom
    1119              :                END IF
    1120              :             END IF
    1121              : 
    1122      8443742 :             IF (PRESENT(nshell)) THEN
    1123        28810 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1124         7425 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, nshell_sum=n)
    1125         7425 :                   nshell = nshell + n*qs_kind_set(ikind)%natom
    1126        21385 :                ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
    1127          960 :                   CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, lmax=n)
    1128          960 :                   nshell = nshell + (n + 1)*qs_kind_set(ikind)%natom
    1129              :                END IF
    1130              :             END IF
    1131              : 
    1132      8443742 :             IF (PRESENT(nelectron)) THEN
    1133       209260 :                IF (ASSOCIATED(qs_kind%all_potential)) THEN
    1134              :                   CALL get_potential(potential=qs_kind%all_potential, &
    1135        18960 :                                      zeff=zeff, zeff_correction=zeff_correction)
    1136       190300 :                ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
    1137              :                   CALL get_potential(potential=qs_kind%gth_potential, &
    1138       188766 :                                      zeff=zeff, zeff_correction=zeff_correction)
    1139         1534 :                ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
    1140              :                   CALL get_potential(potential=qs_kind%sgp_potential, &
    1141          640 :                                      zeff=zeff, zeff_correction=zeff_correction)
    1142              :                ELSE
    1143          894 :                   zeff = 0.0_dp
    1144          894 :                   zeff_correction = 0.0_dp
    1145              :                END IF
    1146       209260 :                nelectron = nelectron + qs_kind_set(ikind)%natom*NINT(zeff - zeff_correction)
    1147              :             END IF
    1148              : 
    1149      8443742 :             IF (PRESENT(basis_rcut)) THEN
    1150          234 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1151            0 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, kind_radius=brcut)
    1152            0 :                   basis_rcut = MAX(basis_rcut, brcut)
    1153          234 :                ELSE IF (ASSOCIATED(qs_kind%dftb_parameter)) THEN
    1154          234 :                   CALL get_dftb_atom_param(dftb_parameter=dftb_parameter, cutoff=brcut)
    1155          234 :                   basis_rcut = MAX(basis_rcut, brcut)
    1156              :                END IF
    1157              :             END IF
    1158              : 
    1159      8443742 :             IF (PRESENT(total_zeff_corr)) THEN
    1160        14319 :                IF (ASSOCIATED(qs_kind%all_potential)) THEN
    1161              :                   CALL get_potential(potential=qs_kind%all_potential, &
    1162         5982 :                                      zeff=zeff, zeff_correction=zeff_correction)
    1163         8337 :                ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
    1164              :                   CALL get_potential(potential=qs_kind%gth_potential, &
    1165         8145 :                                      zeff=zeff, zeff_correction=zeff_correction)
    1166          192 :                ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
    1167              :                   CALL get_potential(potential=qs_kind%sgp_potential, &
    1168           40 :                                      zeff=zeff, zeff_correction=zeff_correction)
    1169              :                ELSE
    1170          152 :                   zeff = 0.0_dp
    1171          152 :                   zeff_correction = 0.0_dp
    1172              :                END IF
    1173        14319 :                total_zeff_corr = total_zeff_corr + qs_kind_set(ikind)%natom*zeff_correction
    1174              :             END IF
    1175              : 
    1176      8443742 :             IF (PRESENT(all_potential_present)) THEN
    1177        67783 :                IF (ASSOCIATED(all_potential)) THEN
    1178        41028 :                   all_potential_present = .TRUE.
    1179              :                END IF
    1180              :             END IF
    1181              : 
    1182      8443742 :             IF (PRESENT(tnadd_potential_present)) THEN
    1183            0 :                IF (ASSOCIATED(tnadd_potential)) THEN
    1184            0 :                   tnadd_potential_present = .TRUE.
    1185              :                END IF
    1186              :             END IF
    1187              : 
    1188      8443742 :             IF (PRESENT(gth_potential_present)) THEN
    1189        53412 :                IF (ASSOCIATED(gth_potential)) THEN
    1190        18172 :                   gth_potential_present = .TRUE.
    1191              :                END IF
    1192              :             END IF
    1193              : 
    1194      8443742 :             IF (PRESENT(sgp_potential_present)) THEN
    1195        53417 :                IF (ASSOCIATED(sgp_potential)) THEN
    1196           60 :                   sgp_potential_present = .TRUE.
    1197              :                END IF
    1198              :             END IF
    1199              : 
    1200      8443742 :             IF (PRESENT(paw_atom_present)) THEN
    1201        52742 :                IF (paw_atom) THEN
    1202         2956 :                   paw_atom_present = .TRUE.
    1203              :                END IF
    1204              :             END IF
    1205              : 
    1206      8443742 :             IF (PRESENT(dft_plus_u_atom_present)) THEN
    1207        14319 :                IF (dft_plus_u_atom) THEN
    1208           32 :                   dft_plus_u_atom_present = .TRUE.
    1209              :                END IF
    1210              :             END IF
    1211              : 
    1212      8443742 :             IF (PRESENT(max_ngrid_rad)) THEN
    1213            0 :                max_ngrid_rad = MAX(max_ngrid_rad, ngrid_rad)
    1214              :             END IF
    1215              : 
    1216      8443742 :             IF (PRESENT(max_sph_harm)) THEN
    1217            0 :                max_sph_harm = MAX(max_sph_harm, max_s_harm)
    1218              :             END IF
    1219              : 
    1220      8443742 :             IF (PRESENT(maxg_iso_not0)) THEN
    1221        31366 :                maxg_iso_not0 = MAX(maxg_iso_not0, max_iso_not0)
    1222              :             END IF
    1223              : 
    1224      8443742 :             IF (PRESENT(lmax_rho0)) THEN
    1225            0 :                lmax_rho0 = MAX(lmax_rho0, lmax_rho0_kind)
    1226              :             END IF
    1227              : 
    1228     20906432 :             IF (PRESENT(npgf_seg)) THEN
    1229           10 :                IF (ASSOCIATED(tmp_basis_set)) THEN
    1230           10 :                   CALL get_gto_basis_set(gto_basis_set=tmp_basis_set, npgf_seg_sum=n)
    1231           10 :                   npgf_seg = npgf_seg + n*qs_kind_set(ikind)%natom
    1232              :                END IF
    1233              :             END IF
    1234              : 
    1235              :          END DO
    1236              :       ELSE
    1237            0 :          CPABORT("The pointer qs_kind_set is not associated")
    1238              :       END IF
    1239              : 
    1240      4018948 :    END SUBROUTINE get_qs_kind_set
    1241              : 
    1242              : ! **************************************************************************************************
    1243              : !> \brief Initialise an atomic kind data set.
    1244              : !> \param qs_kind ...
    1245              : !> \author Creation (11.01.2002,MK)
    1246              : !>                20.09.2002 adapted for pol/kg use, gtb
    1247              : ! **************************************************************************************************
    1248        14319 :    SUBROUTINE init_qs_kind(qs_kind)
    1249              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
    1250              : 
    1251              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'init_qs_kind'
    1252              : 
    1253              :       CHARACTER(LEN=default_string_length)               :: basis_type
    1254              :       INTEGER                                            :: handle, i
    1255              :       TYPE(gto_basis_set_type), POINTER                  :: tmp_basis_set
    1256              : 
    1257        14319 :       CALL timeset(routineN, handle)
    1258              : 
    1259        14319 :       CPASSERT(ASSOCIATED(qs_kind))
    1260              : 
    1261        14319 :       IF (ASSOCIATED(qs_kind%gth_potential)) THEN
    1262         8145 :          CALL init_potential(qs_kind%gth_potential)
    1263         6174 :       ELSEIF (ASSOCIATED(qs_kind%sgp_potential)) THEN
    1264           40 :          CALL init_potential(qs_kind%sgp_potential)
    1265              :       END IF
    1266              : 
    1267       300699 :       DO i = 1, SIZE(qs_kind%basis_sets, 1)
    1268       286380 :          NULLIFY (tmp_basis_set)
    1269              :          CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
    1270       286380 :                                        inumbas=i, basis_type=basis_type)
    1271       286380 :          IF (basis_type == "") CYCLE
    1272        30539 :          IF (basis_type == "AUX") THEN
    1273            0 :             IF (tmp_basis_set%norm_type < 0) tmp_basis_set%norm_type = 1
    1274            0 :             CALL init_aux_basis_set(tmp_basis_set)
    1275              :          ELSE
    1276        16220 :             IF (tmp_basis_set%norm_type < 0) tmp_basis_set%norm_type = 2
    1277        16220 :             CALL init_orb_basis_set(tmp_basis_set)
    1278              :          END IF
    1279              :       END DO
    1280              : 
    1281        14319 :       CALL timestop(handle)
    1282              : 
    1283        14319 :    END SUBROUTINE init_qs_kind
    1284              : 
    1285              : ! **************************************************************************************************
    1286              : !> \brief Initialise an atomic kind set data set.
    1287              : !> \param qs_kind_set ...
    1288              : !> \author - Creation (17.01.2002,MK)
    1289              : !>      - 20.09.2002 para_env passed (gt)
    1290              : ! **************************************************************************************************
    1291         7404 :    SUBROUTINE init_qs_kind_set(qs_kind_set)
    1292              : 
    1293              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1294              : 
    1295              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'init_qs_kind_set'
    1296              : 
    1297              :       INTEGER                                            :: handle, ikind
    1298              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
    1299              : 
    1300         7404 :       CALL timeset(routineN, handle)
    1301              : 
    1302         7404 :       IF (.NOT. ASSOCIATED(qs_kind_set)) THEN
    1303            0 :          CPABORT("init_qs_kind_set: The pointer qs_kind_set is not associated")
    1304              :       END IF
    1305              : 
    1306        21723 :       DO ikind = 1, SIZE(qs_kind_set)
    1307        14319 :          qs_kind => qs_kind_set(ikind)
    1308        21723 :          CALL init_qs_kind(qs_kind)
    1309              :       END DO
    1310              : 
    1311         7404 :       CALL timestop(handle)
    1312              : 
    1313         7404 :    END SUBROUTINE init_qs_kind_set
    1314              : 
    1315              : ! **************************************************************************************************
    1316              : !> \brief ...
    1317              : !> \param qs_kind_set ...
    1318              : !> \param qs_control ...
    1319              : !> \param force_env_section ...
    1320              : !> \param modify_qs_control  whether the qs_control should be modified
    1321              : ! **************************************************************************************************
    1322         1030 :    SUBROUTINE init_gapw_basis_set(qs_kind_set, qs_control, force_env_section, modify_qs_control)
    1323              : 
    1324              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1325              :       TYPE(qs_control_type), POINTER                     :: qs_control
    1326              :       TYPE(section_vals_type), POINTER                   :: force_env_section
    1327              :       LOGICAL, OPTIONAL                                  :: modify_qs_control
    1328              : 
    1329              :       CHARACTER(LEN=default_string_length)               :: bsname
    1330              :       INTEGER                                            :: bas1c, ikind, ilevel, nkind
    1331              :       LOGICAL                                            :: gpw, my_mod_control, paw_atom
    1332              :       REAL(dp)                                           :: max_rad_local_type, rc
    1333              :       TYPE(gto_basis_set_type), POINTER                  :: basis_1c, orb_basis, soft_basis
    1334              :       TYPE(paw_proj_set_type), POINTER                   :: paw_proj
    1335              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
    1336              : 
    1337         1030 :       my_mod_control = .TRUE.
    1338         1030 :       IF (PRESENT(modify_qs_control)) THEN
    1339           92 :          my_mod_control = modify_qs_control
    1340              :       END IF
    1341              : 
    1342         1030 :       IF (ASSOCIATED(qs_kind_set)) THEN
    1343              : 
    1344         1030 :          IF (my_mod_control) qs_control%gapw_control%non_paw_atoms = .FALSE.
    1345         1030 :          nkind = SIZE(qs_kind_set)
    1346              : 
    1347         3008 :          DO ikind = 1, nkind
    1348              : 
    1349         1978 :             qs_kind => qs_kind_set(ikind)
    1350              : 
    1351         1978 :             CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis)
    1352              :             CALL get_qs_kind(qs_kind=qs_kind, hard_radius=rc, &
    1353         1978 :                              max_rad_local=max_rad_local_type, gpw_type_forced=gpw)
    1354              : 
    1355         1978 :             NULLIFY (soft_basis)
    1356         1978 :             CALL allocate_gto_basis_set(soft_basis)
    1357              :             CALL create_soft_basis(orb_basis, soft_basis, &
    1358              :                                    qs_control%gapw_control%eps_fit, rc, paw_atom, &
    1359         1978 :                                    qs_control%gapw_control%force_paw, gpw)
    1360         1978 :             CALL add_basis_set_to_container(qs_kind%basis_sets, soft_basis, "ORB_SOFT")
    1361         1978 :             CALL set_qs_kind(qs_kind=qs_kind, paw_atom=paw_atom)
    1362              : 
    1363         1978 :             bas1c = qs_control%gapw_control%basis_1c
    1364         1978 :             NULLIFY (basis_1c)
    1365         1936 :             SELECT CASE (bas1c)
    1366              :             CASE (gapw_1c_orb)
    1367         1936 :                ilevel = 0
    1368              :             CASE (gapw_1c_small)
    1369           26 :                ilevel = 1
    1370              :             CASE (gapw_1c_medium)
    1371            4 :                ilevel = 2
    1372              :             CASE (gapw_1c_large)
    1373            8 :                ilevel = 3
    1374              :             CASE (gapw_1c_very_large)
    1375            4 :                ilevel = 4
    1376              :             CASE DEFAULT
    1377         1978 :                CPABORT("basis_1c type")
    1378              :             END SELECT
    1379         1978 :             CALL remove_basis_from_container(qs_kind%basis_sets, basis_type="GAPW_1C")
    1380         1978 :             CALL create_1c_basis(orb_basis, soft_basis, basis_1c, ilevel)
    1381         1978 :             CALL get_gto_basis_set(gto_basis_set=orb_basis, name=bsname)
    1382         1978 :             basis_1c%name = TRIM(bsname)//"_1c"
    1383         1978 :             CALL add_basis_set_to_container(qs_kind%basis_sets, basis_1c, "GAPW_1C")
    1384         1978 :             IF (paw_atom) THEN
    1385         1656 :                CALL allocate_paw_proj_set(qs_kind%paw_proj_set)
    1386         1656 :                CALL get_qs_kind(qs_kind=qs_kind, paw_proj_set=paw_proj)
    1387              :                CALL projectors(paw_proj, basis_1c, orb_basis, rc, qs_control, &
    1388         1656 :                                max_rad_local_type, force_env_section)
    1389              :             ELSE
    1390          322 :                IF (my_mod_control) qs_control%gapw_control%non_paw_atoms = .TRUE.
    1391              :             END IF
    1392              : 
    1393              :             ! grid_atom and harmonics are allocated even if NOT PAW_ATOM
    1394         1978 :             NULLIFY (qs_kind%grid_atom, qs_kind%harmonics)
    1395         1978 :             CALL allocate_grid_atom(qs_kind%grid_atom)
    1396         6964 :             CALL allocate_harmonics_atom(qs_kind%harmonics)
    1397              : 
    1398              :          END DO
    1399              : 
    1400         1030 :          IF (my_mod_control) THEN
    1401          938 :             IF (qs_control%gapw_control%non_paw_atoms) THEN
    1402          156 :                qs_control%gapw_control%nopaw_as_gpw = .TRUE.
    1403              :             ELSE
    1404          782 :                qs_control%gapw_control%nopaw_as_gpw = .FALSE.
    1405              :             END IF
    1406              :          END IF
    1407              :       ELSE
    1408            0 :          CPABORT("The pointer qs_kind_set is not associated")
    1409              :       END IF
    1410              : 
    1411         1030 :    END SUBROUTINE init_gapw_basis_set
    1412              : ! **************************************************************************************************
    1413              : !> \brief ...
    1414              : !> \param qs_kind_set ...
    1415              : ! **************************************************************************************************
    1416         1030 :    SUBROUTINE init_gapw_nlcc(qs_kind_set)
    1417              : 
    1418              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    1419              : 
    1420              :       INTEGER                                            :: i, ic, ikind, n_nlcc, nc, nexp_nlcc, &
    1421              :                                                             nkind, nr
    1422         1030 :       INTEGER, DIMENSION(:), POINTER                     :: nct_nlcc
    1423              :       LOGICAL                                            :: nlcc, nlcc_type, paw_atom
    1424              :       REAL(dp)                                           :: alpha, coa, cval
    1425         1030 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: a_nlcc, alpha_nlcc, c_nlcc, fe, rc, rr
    1426         1030 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: cval_nlcc, den
    1427              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
    1428              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
    1429              :       TYPE(sgp_potential_type), POINTER                  :: sgp_potential
    1430              : 
    1431         1030 :       IF (ASSOCIATED(qs_kind_set)) THEN
    1432         1030 :          nlcc = has_nlcc(qs_kind_set)
    1433         1030 :          IF (nlcc) THEN
    1434            2 :             nkind = SIZE(qs_kind_set)
    1435            4 :             DO ikind = 1, nkind
    1436            2 :                qs_kind => qs_kind_set(ikind)
    1437            2 :                CALL get_qs_kind(qs_kind, paw_atom=paw_atom)
    1438            4 :                IF (paw_atom) THEN
    1439            2 :                   CALL get_qs_kind(qs_kind, gth_potential=gth_potential)
    1440            2 :                   CALL get_qs_kind(qs_kind, sgp_potential=sgp_potential)
    1441            2 :                   IF (ASSOCIATED(gth_potential)) THEN
    1442              :                      CALL get_potential(potential=gth_potential, nlcc_present=nlcc_type, &
    1443            2 :                                         nexp_nlcc=nexp_nlcc, alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
    1444            2 :                      IF (nlcc_type) THEN
    1445            2 :                         nr = qs_kind%grid_atom%nr
    1446            2 :                         rr => qs_kind%grid_atom%rad
    1447           12 :                         ALLOCATE (qs_kind%nlcc_pot(nr, 2), rc(nr), fe(nr))
    1448            6 :                         den => qs_kind%nlcc_pot
    1449          206 :                         den = 0.0_dp
    1450            4 :                         DO i = 1, nexp_nlcc
    1451            2 :                            alpha = alpha_nlcc(i)
    1452          202 :                            rc(:) = rr(:)/alpha
    1453          202 :                            fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
    1454            2 :                            nc = nct_nlcc(i)
    1455            8 :                            DO ic = 1, nc
    1456            4 :                               cval = cval_nlcc(ic, i)
    1457            4 :                               coa = cval/alpha
    1458          404 :                               den(:, 1) = den(:, 1) + fe(:)*rc**(2*ic - 2)*cval
    1459          404 :                               den(:, 2) = den(:, 2) - fe(:)*rc**(2*ic - 1)*coa
    1460            6 :                               IF (ic > 1) THEN
    1461          202 :                                  den(:, 2) = den(:, 2) + REAL(2*ic - 2, dp)*fe(:)*rc**(2*ic - 3)*coa
    1462              :                               END IF
    1463              :                            END DO
    1464              :                         END DO
    1465            2 :                         DEALLOCATE (rc, fe)
    1466              :                      END IF
    1467            0 :                   ELSE IF (ASSOCIATED(sgp_potential)) THEN
    1468              :                      CALL get_potential(potential=sgp_potential, has_nlcc=nlcc_type, &
    1469            0 :                                         n_nlcc=n_nlcc, a_nlcc=a_nlcc, c_nlcc=c_nlcc)
    1470            0 :                      IF (nlcc_type) THEN
    1471            0 :                         nr = qs_kind%grid_atom%nr
    1472            0 :                         rr => qs_kind%grid_atom%rad
    1473            0 :                         ALLOCATE (qs_kind%nlcc_pot(nr, 2), rc(nr), fe(nr))
    1474            0 :                         den => qs_kind%nlcc_pot
    1475            0 :                         den = 0.0_dp
    1476            0 :                         DO i = 1, n_nlcc
    1477            0 :                            alpha = a_nlcc(i)
    1478            0 :                            fe(:) = EXP(-alpha*rr(:)*rr(:))
    1479            0 :                            cval = c_nlcc(i)
    1480            0 :                            den(:, 1) = den(:, 1) + cval*fe(:)
    1481            0 :                            den(:, 2) = den(:, 2) - 2.0_dp*alpha*cval*rr(:)*fe(:)
    1482              :                         END DO
    1483            0 :                         DEALLOCATE (rc, fe)
    1484              :                      END IF
    1485              :                   ELSE
    1486              :                      ! skip
    1487              :                   END IF
    1488              :                END IF
    1489              :             END DO
    1490              :          END IF
    1491              :       ELSE
    1492            0 :          CPABORT("The pointer qs_kind_set is not associated")
    1493              :       END IF
    1494              : 
    1495         1030 :    END SUBROUTINE init_gapw_nlcc
    1496              : 
    1497              : ! **************************************************************************************************
    1498              : !> \brief Read an atomic kind data set from the input file.
    1499              : !> \param qs_kind ...
    1500              : !> \param kind_section ...
    1501              : !> \param para_env ...
    1502              : !> \param force_env_section ...
    1503              : !> \param no_fail ...
    1504              : !> \param method_id ...
    1505              : !> \param silent ...
    1506              : !> \par History
    1507              : !>      - Creation (09.02.2002,MK)
    1508              : !>      - 20.09.2002,gt: adapted for POL/KG use (elp_potential)
    1509              : !>      - 05.03.2010: split elp_potential into fist_potential and kg_potential
    1510              : ! **************************************************************************************************
    1511        14397 :    SUBROUTINE read_qs_kind(qs_kind, kind_section, para_env, force_env_section, &
    1512              :                            no_fail, method_id, silent)
    1513              : 
    1514              :       TYPE(qs_kind_type), INTENT(INOUT)                  :: qs_kind
    1515              :       TYPE(section_vals_type), POINTER                   :: kind_section
    1516              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    1517              :       TYPE(section_vals_type), POINTER                   :: force_env_section
    1518              :       LOGICAL, INTENT(IN)                                :: no_fail
    1519              :       INTEGER, INTENT(IN)                                :: method_id
    1520              :       LOGICAL, INTENT(IN)                                :: silent
    1521              : 
    1522              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'read_qs_kind'
    1523              :       INTEGER, PARAMETER                                 :: maxbas = 20
    1524              : 
    1525              :       CHARACTER(LEN=2)                                   :: element_symbol
    1526              :       CHARACTER(len=default_path_length)                 :: kg_potential_fn_kind, &
    1527              :                                                             potential_file_name, potential_fn_kind
    1528              :       CHARACTER(LEN=default_string_length)               :: akind_name, basis_type, keyword, &
    1529              :                                                             kgpot_name, kgpot_type, &
    1530              :                                                             potential_name, potential_type, tmp
    1531              :       CHARACTER(LEN=default_string_length), DIMENSION(4) :: description
    1532              :       CHARACTER(LEN=default_string_length), &
    1533        14397 :          DIMENSION(:), POINTER                           :: tmpstringlist
    1534              :       CHARACTER(LEN=default_string_length), &
    1535              :          DIMENSION(maxbas)                               :: basis_set_form, basis_set_name, &
    1536              :                                                             basis_set_type
    1537              :       INTEGER :: handle, i, i_rep, iounit, ipaodesc, ipaopot, ipos, j, jj, k_rep, l, m, n_rep, &
    1538              :          nb_rep, nexp, ngauss, nlcc, nloc, nnl, norbitals, npaodesc, npaopot, nppnl, nspin, nu, z
    1539        28794 :       INTEGER, DIMENSION(:), POINTER                     :: add_el, elec_conf, orbitals
    1540              :       LOGICAL :: check, ecp_semi_local, explicit, explicit_basis, explicit_J, explicit_kgpot, &
    1541              :          explicit_potential, explicit_U, explicit_u_m_j, nobasis, section_enabled, &
    1542              :          subsection_enabled, update_input
    1543              :       REAL(KIND=dp)                                      :: alpha, ccore, r, rc, zeff_correction
    1544              :       REAL(KIND=dp), DIMENSION(6)                        :: error
    1545        28794 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: a_nl, aloc, anlcc, cloc, cnlcc, nelec
    1546        14397 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: h_nl
    1547        14397 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: c_nl
    1548              :       TYPE(atom_ecppot_type)                             :: ecppot
    1549              :       TYPE(atom_sgp_potential_type)                      :: sgppot
    1550      1511685 :       TYPE(atom_upfpot_type)                             :: upfpot
    1551              :       TYPE(cp_logger_type), POINTER                      :: logger
    1552              :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set, sup_basis_set, &
    1553              :                                                             tmp_basis_set
    1554              :       TYPE(section_vals_type), POINTER :: basis_section, bs_section, dft_plus_u_section, &
    1555              :          dft_section, enforce_occupation_section, kgpot_section, pao_desc_section, &
    1556              :          pao_pot_section, potential_section, spin_section
    1557              :       TYPE(sto_basis_set_type), POINTER                  :: sto_basis_set
    1558              : 
    1559        14397 :       CALL timeset(routineN, handle)
    1560              : 
    1561        14397 :       NULLIFY (logger)
    1562        14397 :       logger => cp_get_default_logger()
    1563        14397 :       iounit = cp_logger_get_default_io_unit(logger)
    1564              : 
    1565        14397 :       NULLIFY (elec_conf)
    1566              : 
    1567        14397 :       update_input = .TRUE.
    1568       302337 :       basis_set_name(:) = ""
    1569       302337 :       basis_set_type(:) = ""
    1570       302337 :       basis_set_form(:) = ""
    1571        14397 :       potential_name = ""
    1572        14397 :       potential_type = ""
    1573        14397 :       kgpot_name = ""
    1574        14397 :       kgpot_type = ""
    1575        14397 :       z = -1
    1576        14397 :       zeff_correction = 0.0_dp
    1577        14397 :       explicit = .FALSE.
    1578        14397 :       explicit_basis = .FALSE.
    1579        14397 :       explicit_J = .FALSE.
    1580        14397 :       explicit_kgpot = .FALSE.
    1581        14397 :       explicit_potential = .FALSE.
    1582        14397 :       explicit_U = .FALSE.
    1583        14397 :       explicit_u_m_j = .FALSE.
    1584              : 
    1585        14397 :       dft_section => section_vals_get_subs_vals(force_env_section, "DFT")
    1586        14397 :       CALL section_vals_get(kind_section, n_repetition=n_rep)
    1587        14397 :       k_rep = -1
    1588        14397 :       akind_name = qs_kind%name
    1589        14397 :       CALL uppercase(akind_name)
    1590              :       ! First we use the atom_name to find out the proper KIND section
    1591        20916 :       DO i_rep = 1, n_rep
    1592              :          CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
    1593        16028 :                                    c_val=keyword, i_rep_section=i_rep)
    1594        16028 :          CALL uppercase(keyword)
    1595        20916 :          IF (keyword == akind_name) THEN
    1596         9509 :             k_rep = i_rep
    1597         9509 :             EXIT
    1598              :          END IF
    1599              :       END DO
    1600              :       ! The search for the KIND section failed.. check for a QM/MM link atom
    1601        14397 :       IF (k_rep < 1) THEN
    1602         4888 :          ipos = INDEX(qs_kind%name, "_")
    1603         4888 :          IF (((ipos == 2) .OR. (ipos == 3)) .AND. (INDEX(qs_kind%name, "_ghost") == 0)) THEN
    1604              :             ! If the atm_name could not match any KIND section it maybe be a QM/MM link atom.
    1605              :             ! ghost atoms will be treated differently.
    1606           64 :             akind_name = qs_kind%name(1:ipos - 1)
    1607           64 :             CALL uppercase(akind_name)
    1608           64 :             DO i_rep = 1, n_rep
    1609              :                CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
    1610           52 :                                          c_val=keyword, i_rep_section=i_rep)
    1611           52 :                CALL uppercase(keyword)
    1612           64 :                IF (keyword == akind_name) THEN
    1613           52 :                   k_rep = i_rep
    1614           52 :                   EXIT
    1615              :                END IF
    1616              :             END DO
    1617              :          END IF
    1618              :       END IF
    1619              :       ! The search for the KIND section failed.. check element_symbol
    1620        14397 :       IF (k_rep < 1) THEN
    1621              :          ! If it's not a link atom let's check for the element and map
    1622              :          ! the KIND section to the element.
    1623         4836 :          element_symbol = qs_kind%element_symbol(1:2)
    1624         4836 :          CALL uppercase(element_symbol)
    1625         4928 :          DO i_rep = 1, n_rep
    1626              :             CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
    1627          112 :                                       c_val=keyword, i_rep_section=i_rep)
    1628          112 :             CALL uppercase(keyword)
    1629         4928 :             IF (keyword == element_symbol) THEN
    1630           20 :                k_rep = i_rep
    1631           20 :                EXIT
    1632              :             END IF
    1633              :          END DO
    1634              :       END IF
    1635              :       ! In case it should not really match any possible KIND section
    1636              :       ! let's look if a default one is defined..
    1637        14397 :       IF (k_rep < 1) THEN
    1638         4832 :          DO i_rep = 1, n_rep
    1639              :             CALL section_vals_val_get(kind_section, "_SECTION_PARAMETERS_", &
    1640           72 :                                       c_val=keyword, i_rep_section=i_rep)
    1641           72 :             CALL uppercase(keyword)
    1642         4832 :             IF (keyword == "DEFAULT") THEN
    1643           56 :                update_input = .FALSE.
    1644           56 :                k_rep = i_rep
    1645           56 :                EXIT
    1646              :             END IF
    1647              :          END DO
    1648              :       END IF
    1649        14397 :       IF (k_rep < 0 .AND. (.NOT. no_fail)) THEN
    1650              :          CALL cp_abort(__LOCATION__, &
    1651              :                        "No &KIND section was possible to associate to the atomic kind <"// &
    1652              :                        TRIM(akind_name)//">. The KIND section were also scanned for the"// &
    1653              :                        " corresponding element <"//TRIM(qs_kind%element_symbol)//">"// &
    1654            0 :                        " and for the DEFAULT section but no match was found. Check your input file!")
    1655              :       END IF
    1656              :       ! Retrieve information on element
    1657        14397 :       CALL get_ptable_info(qs_kind%element_symbol, ielement=z)
    1658              : 
    1659              :       ! Normal parsing of the KIND section
    1660        14397 :       IF (k_rep > 0) THEN
    1661              :          ! new style basis set input
    1662              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1663              :                                    keyword_name="BASIS_SET", &
    1664              :                                    explicit=explicit, &
    1665         9637 :                                    n_rep_val=nb_rep)
    1666         9637 :          IF (.NOT. explicit) nb_rep = 0
    1667         9637 :          CPASSERT(nb_rep <= maxbas)
    1668        21085 :          DO i = 1, nb_rep
    1669              :             CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1670        11448 :                                       keyword_name="BASIS_SET", i_rep_val=i, c_vals=tmpstringlist)
    1671        11448 :             IF (SIZE(tmpstringlist) == 1) THEN
    1672              :                ! default is orbital type and GTO
    1673         8635 :                basis_set_type(i) = "ORB"
    1674         8635 :                basis_set_form(i) = "GTO"
    1675         8635 :                basis_set_name(i) = tmpstringlist(1)
    1676         2813 :             ELSEIF (SIZE(tmpstringlist) == 2) THEN
    1677              :                ! default is GTO
    1678         2809 :                basis_set_type(i) = tmpstringlist(1)
    1679         2809 :                basis_set_form(i) = "GTO"
    1680         2809 :                basis_set_name(i) = tmpstringlist(2)
    1681            4 :             ELSEIF (SIZE(tmpstringlist) == 3) THEN
    1682            4 :                basis_set_type(i) = tmpstringlist(1)
    1683            4 :                basis_set_form(i) = tmpstringlist(2)
    1684            4 :                basis_set_name(i) = tmpstringlist(3)
    1685              :             ELSE
    1686              :                CALL cp_abort(__LOCATION__, &
    1687            0 :                              "invalid number of BASIS_SET keyword parameters: BASIS_SET [<TYPE>] [<FORM>] <NAME>")
    1688              :             END IF
    1689              :             ! check that we have a valid basis set form
    1690        21085 :             IF (basis_set_form(i) /= "GTO" .AND. basis_set_form(i) /= "STO") THEN
    1691            0 :                CPABORT("invalid BASIS_SET FORM parameter")
    1692              :             END IF
    1693              :          END DO
    1694              : 
    1695              :          ! parse PAO keywords
    1696              :          CALL section_vals_val_get(kind_section, keyword_name="PAO_BASIS_SIZE", i_rep_section=k_rep, &
    1697         9637 :                                    i_val=qs_kind%pao_basis_size)
    1698              :          CALL section_vals_val_get(kind_section, keyword_name="PAO_MODEL_FILE", i_rep_section=k_rep, &
    1699         9637 :                                    explicit=explicit)
    1700         9637 :          IF (explicit) THEN
    1701              :             CALL section_vals_val_get(kind_section, keyword_name="PAO_MODEL_FILE", i_rep_section=k_rep, &
    1702            8 :                                       c_val=qs_kind%pao_model_file)
    1703              :          END IF
    1704              : 
    1705              :          ! parse PAO_POTENTIAL sections
    1706         9637 :          pao_pot_section => section_vals_get_subs_vals(kind_section, "PAO_POTENTIAL", i_rep_section=k_rep)
    1707         9637 :          CALL section_vals_get(pao_pot_section, n_repetition=npaopot)
    1708        19396 :          ALLOCATE (qs_kind%pao_potentials(npaopot))
    1709         9699 :          DO ipaopot = 1, npaopot
    1710              :             CALL section_vals_val_get(pao_pot_section, keyword_name="MAXL", i_rep_section=ipaopot, &
    1711           62 :                                       i_val=qs_kind%pao_potentials(ipaopot)%maxl)
    1712              :             CALL section_vals_val_get(pao_pot_section, keyword_name="MAX_PROJECTOR", i_rep_section=ipaopot, &
    1713           62 :                                       i_val=qs_kind%pao_potentials(ipaopot)%max_projector)
    1714              :             CALL section_vals_val_get(pao_pot_section, keyword_name="BETA", i_rep_section=ipaopot, &
    1715           62 :                                       r_val=qs_kind%pao_potentials(ipaopot)%beta)
    1716              :             CALL section_vals_val_get(pao_pot_section, keyword_name="WEIGHT", i_rep_section=ipaopot, &
    1717         9699 :                                       r_val=qs_kind%pao_potentials(ipaopot)%weight)
    1718              :          END DO
    1719              : 
    1720              :          ! parse PAO_DESCRIPTOR sections
    1721         9637 :          pao_desc_section => section_vals_get_subs_vals(kind_section, "PAO_DESCRIPTOR", i_rep_section=k_rep)
    1722         9637 :          CALL section_vals_get(pao_desc_section, n_repetition=npaodesc)
    1723        19304 :          ALLOCATE (qs_kind%pao_descriptors(npaodesc))
    1724         9655 :          DO ipaodesc = 1, npaodesc
    1725              :             CALL section_vals_val_get(pao_desc_section, keyword_name="BETA", i_rep_section=ipaodesc, &
    1726           18 :                                       r_val=qs_kind%pao_descriptors(ipaodesc)%beta)
    1727              :             CALL section_vals_val_get(pao_desc_section, keyword_name="SCREENING", i_rep_section=ipaodesc, &
    1728           18 :                                       r_val=qs_kind%pao_descriptors(ipaodesc)%screening)
    1729              :             CALL section_vals_val_get(pao_desc_section, keyword_name="WEIGHT", i_rep_section=ipaodesc, &
    1730         9655 :                                       r_val=qs_kind%pao_descriptors(ipaodesc)%weight)
    1731              :          END DO
    1732              : 
    1733              :          ! parse ELEC_CONF
    1734              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1735         9637 :                                    keyword_name="ELEC_CONF", n_rep_val=i)
    1736         9637 :          IF (i > 0) THEN
    1737              :             CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1738            4 :                                       keyword_name="ELEC_CONF", i_vals=elec_conf)
    1739            4 :             CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
    1740              :          END IF
    1741              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1742         9637 :                                    keyword_name="CORE_CORRECTION", r_val=zeff_correction)
    1743              :          ! parse POTENTIAL
    1744              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1745         9637 :                                    keyword_name="POTENTIAL_FILE_NAME", c_val=potential_fn_kind)
    1746              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1747         9637 :                                    keyword_name="POTENTIAL_TYPE", c_val=potential_type)
    1748              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1749         9637 :                                    explicit=explicit, keyword_name="POTENTIAL", c_vals=tmpstringlist)
    1750         9637 :          IF (explicit) THEN
    1751         9361 :             IF (SIZE(tmpstringlist) == 1) THEN
    1752              :                ! old type of input: start of name defines type
    1753         9301 :                potential_name = tmpstringlist(1)
    1754         9301 :                IF (potential_type == "") THEN
    1755         9301 :                   ipos = INDEX(potential_name, "-")
    1756         9301 :                   IF (ipos > 1) THEN
    1757         8241 :                      potential_type = potential_name(:ipos - 1)
    1758              :                   ELSE
    1759         1060 :                      potential_type = potential_name
    1760              :                   END IF
    1761              :                END IF
    1762           60 :             ELSEIF (SIZE(tmpstringlist) == 2) THEN
    1763           60 :                potential_type = tmpstringlist(1)
    1764           60 :                potential_name = tmpstringlist(2)
    1765              :             ELSE
    1766            0 :                CPABORT("POTENTIAL input list is not correct")
    1767              :             END IF
    1768              :          END IF
    1769         9637 :          CALL uppercase(potential_type)
    1770              : 
    1771              :          ! Parse KG POTENTIAL
    1772              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1773         9637 :                                    keyword_name="KG_POTENTIAL_FILE_NAME", c_val=kg_potential_fn_kind)
    1774              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1775         9637 :                                    keyword_name="KG_POTENTIAL", c_val=kgpot_name)
    1776              : 
    1777              :          ! Semi-local vs. full nonlocal form of ECPs
    1778              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1779         9637 :                                    keyword_name="ECP_SEMI_LOCAL", l_val=ecp_semi_local)
    1780              : 
    1781              :          ! Assign atomic covalent radius
    1782         9637 :          qs_kind%covalent_radius = ptable(z)%covalent_radius*bohr
    1783              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1784         9637 :                                    keyword_name="COVALENT_RADIUS", r_val=r)
    1785         9637 :          IF (r > 0.0_dp) qs_kind%covalent_radius = r
    1786              : 
    1787              :          ! Assign atomic van der Waals radius
    1788         9637 :          qs_kind%vdw_radius = ptable(z)%vdw_radius*bohr
    1789              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1790         9637 :                                    keyword_name="VDW_RADIUS", r_val=r)
    1791         9637 :          IF (r > 0.0_dp) qs_kind%vdw_radius = r
    1792              : 
    1793              :          ! Assign atom dependent defaults, only H special case
    1794              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, n_rep_val=i, &
    1795         9637 :                                    keyword_name="HARD_EXP_RADIUS")
    1796         9637 :          IF (i == 0) THEN
    1797         9583 :             IF (z == 1) THEN
    1798         4172 :                qs_kind%hard_radius = 1.2_dp
    1799              :             ELSE
    1800         5411 :                qs_kind%hard_radius = 0.8_dp*bohr
    1801              :             END IF
    1802              :          ELSE
    1803              :             CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1804           54 :                                       keyword_name="HARD_EXP_RADIUS", r_val=qs_kind%hard_radius)
    1805              :          END IF
    1806              : 
    1807              :          ! assign atom dependent defaults, only H special case
    1808              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, n_rep_val=i, &
    1809         9637 :                                    keyword_name="RHO0_EXP_RADIUS")
    1810         9637 :          IF (i == 0) THEN
    1811         9637 :             qs_kind%hard0_radius = qs_kind%hard_radius
    1812              :          ELSE
    1813              :             CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1814            0 :                                       keyword_name="RHO0_EXP_RADIUS", r_val=qs_kind%hard0_radius)
    1815              :          END IF
    1816         9637 :          IF (qs_kind%hard_radius < qs_kind%hard0_radius) &
    1817            0 :             CPABORT("rc0 should be <= rc")
    1818              : 
    1819              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1820         9637 :                                    keyword_name="MAX_RAD_LOCAL", r_val=qs_kind%max_rad_local)
    1821              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1822         9637 :                                    keyword_name="LEBEDEV_GRID", i_val=qs_kind%ngrid_ang)
    1823         9637 :          IF (qs_kind%ngrid_ang <= 0) &
    1824            0 :             CPABORT("# point lebedev grid < 0")
    1825              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1826         9637 :                                    keyword_name="RADIAL_GRID", i_val=qs_kind%ngrid_rad)
    1827         9637 :          IF (qs_kind%ngrid_rad <= 0) &
    1828            0 :             CPABORT("# point radial grid < 0")
    1829              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1830         9637 :                                    keyword_name="GPW_TYPE", l_val=qs_kind%gpw_type_forced)
    1831              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1832         9637 :                                    keyword_name="GHOST", l_val=qs_kind%ghost)
    1833              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1834         9637 :                                    keyword_name="FLOATING_BASIS_CENTER", l_val=qs_kind%floating)
    1835              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1836         9637 :                                    keyword_name="NO_OPTIMIZE", l_val=qs_kind%no_optimize)
    1837              : 
    1838              :          ! Magnetization
    1839              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1840         9637 :                                    keyword_name="MAGNETIZATION", r_val=qs_kind%magnetization)
    1841              :          ! DFTB3 param
    1842              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1843         9637 :                                    keyword_name="DFTB3_PARAM", r_val=qs_kind%dudq_dftb3)
    1844              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1845         9637 :                                    keyword_name="LMAX_DFTB", i_val=qs_kind%lmax_dftb)
    1846              : 
    1847              :          ! MAOS
    1848              :          CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    1849         9637 :                                    keyword_name="MAO", i_val=qs_kind%mao)
    1850              : 
    1851              :          ! Read the BS subsection of the current atomic kind, if enabled
    1852         9637 :          NULLIFY (bs_section)
    1853              :          bs_section => section_vals_get_subs_vals(kind_section, "BS", &
    1854         9637 :                                                   i_rep_section=k_rep)
    1855              :          section_enabled = .FALSE.
    1856              :          CALL section_vals_val_get(bs_section, "_SECTION_PARAMETERS_", &
    1857         9637 :                                    l_val=section_enabled)
    1858         9637 :          IF (section_enabled) THEN
    1859              :             ! test for conflict with magnetization
    1860           62 :             IF (qs_kind%magnetization /= 0.0_dp) THEN
    1861              :                CALL cp_abort(__LOCATION__, "BS Section is in conflict with non-zero magnetization "// &
    1862            0 :                              "for this atom kind.")
    1863              :             END IF
    1864           62 :             qs_kind%bs_occupation = .TRUE.
    1865              :             !Alpha spin
    1866           62 :             NULLIFY (spin_section)
    1867           62 :             spin_section => section_vals_get_subs_vals(bs_section, "ALPHA")
    1868           62 :             CALL section_vals_get(spin_section, explicit=explicit)
    1869           62 :             IF (explicit) THEN
    1870           62 :                NULLIFY (add_el)
    1871              :                CALL section_vals_val_get(spin_section, &
    1872           62 :                                          keyword_name="NEL", i_vals=add_el)
    1873           62 :                CPASSERT(ASSOCIATED(add_el))
    1874          186 :                ALLOCATE (qs_kind%addel(SIZE(add_el), 2))
    1875          338 :                qs_kind%addel = 0
    1876          138 :                qs_kind%addel(1:SIZE(add_el), 1) = add_el(1:SIZE(add_el))
    1877           62 :                NULLIFY (add_el)
    1878              :                CALL section_vals_val_get(spin_section, &
    1879           62 :                                          keyword_name="L", i_vals=add_el)
    1880           62 :                CPASSERT(ASSOCIATED(add_el))
    1881           62 :                CPASSERT(SIZE(add_el) == SIZE(qs_kind%addel, 1))
    1882          186 :                ALLOCATE (qs_kind%laddel(SIZE(add_el), 2))
    1883          338 :                qs_kind%laddel = 0
    1884          138 :                qs_kind%laddel(1:SIZE(add_el), 1) = add_el(1:SIZE(add_el))
    1885          186 :                ALLOCATE (qs_kind%naddel(SIZE(add_el), 2))
    1886          338 :                qs_kind%naddel = 0
    1887           62 :                NULLIFY (add_el)
    1888              :                CALL section_vals_val_get(spin_section, &
    1889           62 :                                          keyword_name="N", n_rep_val=i)
    1890           62 :                IF (i > 0) THEN
    1891              :                   CALL section_vals_val_get(spin_section, &
    1892           62 :                                             keyword_name="N", i_vals=add_el)
    1893           62 :                   IF (SIZE(add_el) == SIZE(qs_kind%addel, 1)) THEN
    1894          138 :                      qs_kind%naddel(1:SIZE(add_el), 1) = add_el(1:SIZE(add_el))
    1895              :                   END IF
    1896              :                END IF
    1897              :             END IF
    1898              :             ! Beta spin
    1899           62 :             NULLIFY (spin_section)
    1900           62 :             spin_section => section_vals_get_subs_vals(bs_section, "BETA")
    1901           62 :             CALL section_vals_get(spin_section, explicit=explicit)
    1902           62 :             IF (explicit) THEN
    1903           62 :                NULLIFY (add_el)
    1904              :                CALL section_vals_val_get(spin_section, &
    1905           62 :                                          keyword_name="NEL", i_vals=add_el)
    1906           62 :                CPASSERT(SIZE(add_el) == SIZE(qs_kind%addel, 1))
    1907          138 :                qs_kind%addel(1:SIZE(add_el), 2) = add_el(1:SIZE(add_el))
    1908          338 :                qs_kind%addel(:, :) = qs_kind%addel(:, :)
    1909           62 :                NULLIFY (add_el)
    1910              :                CALL section_vals_val_get(spin_section, &
    1911           62 :                                          keyword_name="L", i_vals=add_el)
    1912           62 :                CPASSERT(SIZE(add_el) == SIZE(qs_kind%addel, 1))
    1913          138 :                qs_kind%laddel(1:SIZE(add_el), 2) = add_el(1:SIZE(add_el))
    1914              : 
    1915              :                CALL section_vals_val_get(spin_section, &
    1916           62 :                                          keyword_name="N", n_rep_val=i)
    1917           62 :                IF (i > 0) THEN
    1918           62 :                   NULLIFY (add_el)
    1919              :                   CALL section_vals_val_get(spin_section, &
    1920           62 :                                             keyword_name="N", i_vals=add_el)
    1921           62 :                   IF (SIZE(add_el) == SIZE(qs_kind%addel, 1)) THEN
    1922          138 :                      qs_kind%naddel(1:SIZE(add_el), 2) = add_el(1:SIZE(add_el))
    1923              :                   END IF
    1924              :                END IF
    1925              :             END IF
    1926              :          END IF
    1927              : 
    1928              :          ! Read the DFT+U subsection of the current atomic kind, if enabled
    1929              : 
    1930         9637 :          NULLIFY (dft_plus_u_section)
    1931              :          dft_plus_u_section => section_vals_get_subs_vals(kind_section, &
    1932              :                                                           subsection_name="DFT_PLUS_U", &
    1933         9637 :                                                           i_rep_section=k_rep)
    1934              :          section_enabled = .FALSE.
    1935              :          CALL section_vals_val_get(dft_plus_u_section, &
    1936              :                                    keyword_name="_SECTION_PARAMETERS_", &
    1937         9637 :                                    l_val=section_enabled)
    1938       125281 :          IF (section_enabled) THEN
    1939           32 :             ALLOCATE (qs_kind%dft_plus_u)
    1940              :             NULLIFY (qs_kind%dft_plus_u%nelec)
    1941              :             NULLIFY (qs_kind%dft_plus_u%orbitals)
    1942              :             CALL section_vals_val_get(dft_plus_u_section, &
    1943              :                                       keyword_name="L", &
    1944           32 :                                       i_val=l)
    1945           32 :             qs_kind%dft_plus_u%l = l
    1946              : #if defined(__SIRIUS)
    1947              :             CALL section_vals_val_get(dft_plus_u_section, &
    1948              :                                       keyword_name="N", &
    1949           32 :                                       i_val=nu)
    1950           32 :             qs_kind%dft_plus_u%n = nu
    1951              : 
    1952              :             CALL section_vals_val_get(dft_plus_u_section, &
    1953              :                                       keyword_name="U", &
    1954              :                                       r_val=qs_kind%dft_plus_u%U, &
    1955           32 :                                       explicit=explicit_U)
    1956              : 
    1957              :             CALL section_vals_val_get(dft_plus_u_section, &
    1958              :                                       keyword_name="J", &
    1959              :                                       r_val=qs_kind%dft_plus_u%J, &
    1960           32 :                                       explicit=explicit_J)
    1961              : 
    1962              :             CALL section_vals_val_get(dft_plus_u_section, &
    1963              :                                       keyword_name="alpha", &
    1964           32 :                                       r_val=qs_kind%dft_plus_u%alpha)
    1965              : 
    1966              :             CALL section_vals_val_get(dft_plus_u_section, &
    1967              :                                       keyword_name="beta", &
    1968           32 :                                       r_val=qs_kind%dft_plus_u%beta)
    1969              : 
    1970              :             CALL section_vals_val_get(dft_plus_u_section, &
    1971              :                                       keyword_name="J0", &
    1972           32 :                                       r_val=qs_kind%dft_plus_u%J0)
    1973              : 
    1974              :             CALL section_vals_val_get(dft_plus_u_section, &
    1975              :                                       keyword_name="occupation", &
    1976           32 :                                       r_val=qs_kind%dft_plus_u%occupation)
    1977              : #else
    1978              :             nu = 0
    1979              : #endif
    1980              : 
    1981              :             CALL section_vals_val_get(dft_plus_u_section, &
    1982              :                                       keyword_name="U_MINUS_J", &
    1983              :                                       r_val=qs_kind%dft_plus_u%u_minus_j_target, &
    1984           32 :                                       explicit=explicit_u_m_j)
    1985              : 
    1986           32 :             IF ((explicit_U .OR. explicit_J) .AND. explicit_u_m_j) THEN
    1987            0 :                CPABORT("DFT+U| specifying U or J and U_MINUS_J parameters are mutually exclusive.")
    1988              :             END IF
    1989              : 
    1990              :             CALL section_vals_val_get(dft_plus_u_section, &
    1991              :                                       keyword_name="U_RAMPING", &
    1992           32 :                                       r_val=qs_kind%dft_plus_u%u_ramping)
    1993              :             CALL section_vals_val_get(dft_plus_u_section, &
    1994              :                                       keyword_name="INIT_U_RAMPING_EACH_SCF", &
    1995           32 :                                       l_val=qs_kind%dft_plus_u%init_u_ramping_each_scf)
    1996           32 :             IF (qs_kind%dft_plus_u%u_ramping > 0.0_dp) THEN
    1997            8 :                qs_kind%dft_plus_u%u_minus_j = 0.0_dp
    1998              :             ELSE
    1999           24 :                qs_kind%dft_plus_u%u_minus_j = qs_kind%dft_plus_u%u_minus_j_target
    2000              :             END IF
    2001              :             CALL section_vals_val_get(dft_plus_u_section, &
    2002              :                                       keyword_name="EPS_U_RAMPING", &
    2003           32 :                                       r_val=qs_kind%dft_plus_u%eps_u_ramping)
    2004              : 
    2005           32 :             NULLIFY (enforce_occupation_section)
    2006              :             enforce_occupation_section => section_vals_get_subs_vals(dft_plus_u_section, &
    2007           32 :                                                                      subsection_name="ENFORCE_OCCUPATION")
    2008              :             subsection_enabled = .FALSE.
    2009              :             CALL section_vals_val_get(enforce_occupation_section, &
    2010              :                                       keyword_name="_SECTION_PARAMETERS_", &
    2011           32 :                                       l_val=subsection_enabled)
    2012           32 :             IF (subsection_enabled) THEN
    2013            4 :                NULLIFY (nelec)
    2014              :                CALL section_vals_val_get(enforce_occupation_section, &
    2015              :                                          keyword_name="NELEC", &
    2016            4 :                                          r_vals=nelec)
    2017            4 :                nspin = SIZE(nelec)
    2018           12 :                ALLOCATE (qs_kind%dft_plus_u%nelec(nspin))
    2019            8 :                qs_kind%dft_plus_u%nelec(:) = nelec(:)
    2020            4 :                NULLIFY (orbitals)
    2021              :                CALL section_vals_val_get(enforce_occupation_section, &
    2022              :                                          keyword_name="ORBITALS", &
    2023            4 :                                          i_vals=orbitals)
    2024            4 :                norbitals = SIZE(orbitals)
    2025            4 :                IF (norbitals <= 0 .OR. norbitals > 2*l + 1) &
    2026              :                   CALL cp_abort(__LOCATION__, "DFT+U| Invalid number of ORBITALS specified: "// &
    2027            0 :                                 "1 to 2*L+1 integer numbers are expected")
    2028           12 :                ALLOCATE (qs_kind%dft_plus_u%orbitals(norbitals))
    2029           16 :                qs_kind%dft_plus_u%orbitals(:) = orbitals(:)
    2030            4 :                NULLIFY (orbitals)
    2031           16 :                DO m = 1, norbitals
    2032           12 :                   IF (qs_kind%dft_plus_u%orbitals(m) > l) &
    2033            0 :                      CPABORT("DFT+U| Invalid orbital magnetic quantum number specified: m > l")
    2034           12 :                   IF (qs_kind%dft_plus_u%orbitals(m) < -l) &
    2035            0 :                      CPABORT("DFT+U| Invalid orbital magnetic quantum number specified: m < -l")
    2036           52 :                   DO j = 1, norbitals
    2037           48 :                      IF (j /= m) THEN
    2038           24 :                         IF (qs_kind%dft_plus_u%orbitals(j) == qs_kind%dft_plus_u%orbitals(m)) &
    2039            0 :                            CPABORT("DFT+U| An orbital magnetic quantum number was specified twice")
    2040              :                      END IF
    2041              :                   END DO
    2042              :                END DO
    2043              :                CALL section_vals_val_get(enforce_occupation_section, &
    2044              :                                          keyword_name="EPS_SCF", &
    2045            4 :                                          r_val=qs_kind%dft_plus_u%eps_scf)
    2046              :                CALL section_vals_val_get(enforce_occupation_section, &
    2047              :                                          keyword_name="MAX_SCF", &
    2048            4 :                                          i_val=i)
    2049            4 :                qs_kind%dft_plus_u%max_scf = MAX(-1, i)
    2050              :                CALL section_vals_val_get(enforce_occupation_section, &
    2051              :                                          keyword_name="SMEAR", &
    2052            4 :                                          l_val=qs_kind%dft_plus_u%smear)
    2053              :             END IF ! subsection enabled
    2054              :          END IF ! section enabled
    2055              : 
    2056              :       END IF
    2057              : 
    2058              :       ! Allocate and initialise the orbital basis set data set structure
    2059        14397 :       CALL init_orbital_pointers(5) ! debug the SUN optimizer
    2060              : 
    2061              :       ! BASIS  and POTENTIAL read only when strictly necessary otherwise, even if not used
    2062              :       ! we just print misleading informations
    2063        14397 :       explicit_basis = .FALSE.
    2064        14397 :       IF (k_rep > 0) THEN
    2065              :          basis_section => section_vals_get_subs_vals(kind_section, "BASIS", i_rep_section=k_rep, &
    2066         9637 :                                                      can_return_null=.TRUE.)
    2067         9637 :          CALL section_vals_get(basis_section, explicit=explicit_basis)
    2068              :       END IF
    2069              : 
    2070        14397 :       explicit_potential = .FALSE.
    2071        14397 :       IF (k_rep > 0) THEN
    2072              :          potential_section => section_vals_get_subs_vals(kind_section, "POTENTIAL", &
    2073         9637 :                                                          i_rep_section=k_rep, can_return_null=.TRUE.)
    2074         9637 :          CALL section_vals_get(potential_section, explicit=explicit_potential)
    2075              :       END IF
    2076              : 
    2077        14397 :       explicit_kgpot = .FALSE.
    2078        14397 :       IF (k_rep > 0) THEN
    2079              :          kgpot_section => section_vals_get_subs_vals(kind_section, "KG_POTENTIAL", &
    2080         9637 :                                                      i_rep_section=k_rep, can_return_null=.TRUE.)
    2081         9637 :          CALL section_vals_get(kgpot_section, explicit=explicit_kgpot)
    2082              :       END IF
    2083              : 
    2084        16637 :       SELECT CASE (method_id)
    2085              :       CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pdg, do_method_pm3, do_method_pm6, &
    2086              :             do_method_pm6fm, do_method_mndod, do_method_pnnl)
    2087              :          ! Allocate all_potential
    2088         2240 :          CALL allocate_potential(qs_kind%all_potential)
    2089         2240 :          CALL set_default_all_potential(qs_kind%all_potential, z, zeff_correction)
    2090         2240 :          CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
    2091         2240 :          IF (.NOT. ASSOCIATED(elec_conf)) THEN
    2092         2240 :             CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
    2093         2240 :             CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
    2094              :          END IF
    2095         2240 :          CPASSERT(.NOT. qs_kind%floating)
    2096         2240 :          IF (qs_kind%ghost) THEN
    2097            0 :             CALL get_qs_kind(qs_kind=qs_kind, elec_conf=elec_conf)
    2098            0 :             elec_conf(:) = 0
    2099              :             CALL get_potential(potential=qs_kind%all_potential, &
    2100            0 :                                elec_conf=elec_conf)
    2101            0 :             elec_conf(:) = 0
    2102              :             CALL set_potential(potential=qs_kind%all_potential, &
    2103              :                                zeff=0.0_dp, &
    2104            0 :                                zeff_correction=0.0_dp)
    2105              :          END IF
    2106              : 
    2107              :          ! Basis set (Parameters)
    2108              :          ! Setup proper semiempirical parameters
    2109         2240 :          check = .NOT. ASSOCIATED(qs_kind%se_parameter)
    2110         2240 :          CPASSERT(check)
    2111         2240 :          CALL semi_empirical_create(qs_kind%se_parameter)
    2112              :          ! Check if we allow p-orbitals on H
    2113          438 :          SELECT CASE (z)
    2114              :          CASE (1)
    2115         2240 :             IF (k_rep > 0) THEN
    2116              :                CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
    2117           52 :                                          keyword_name="SE_P_ORBITALS_ON_H", l_val=qs_kind%se_parameter%p_orbitals_on_h)
    2118              :             END IF
    2119              :          CASE DEFAULT
    2120              :             ! No special cases for other elements..
    2121              :          END SELECT
    2122              :          ! Set default parameters
    2123         2240 :          CALL section_vals_val_get(dft_section, "QS%SE%STO_NG", i_val=ngauss)
    2124         2240 :          CALL se_param_set_default(qs_kind%se_parameter, z, method_id)
    2125         2240 :          NULLIFY (tmp_basis_set)
    2126         2240 :          CALL init_se_param(qs_kind%se_parameter, tmp_basis_set, ngauss)
    2127         2240 :          CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, "ORB")
    2128              :          CALL init_potential(qs_kind%all_potential, itype="BARE", &
    2129         2240 :                              zeff=qs_kind%se_parameter%zeff, zeff_correction=zeff_correction)
    2130         2240 :          qs_kind%se_parameter%zeff = qs_kind%se_parameter%zeff - zeff_correction
    2131              : 
    2132         2240 :          check = ((potential_name /= '') .OR. explicit_potential) .AND. .NOT. silent
    2133              :          IF (check) &
    2134              :             CALL cp_warn(__LOCATION__, &
    2135              :                          "Information provided in the input file regarding POTENTIAL for KIND <"// &
    2136           80 :                          TRIM(qs_kind%name)//"> will be ignored!")
    2137              : 
    2138         2240 :          check = ((k_rep > 0) .OR. explicit_basis) .AND. .NOT. silent
    2139              :          IF (check) &
    2140              :             CALL cp_warn(__LOCATION__, &
    2141              :                          "Information provided in the input file regarding BASIS for KIND <"// &
    2142          116 :                          TRIM(qs_kind%name)//"> will be ignored!")
    2143              : 
    2144              :       CASE (do_method_dftb)
    2145              :          ! Allocate all_potential
    2146          480 :          CALL allocate_potential(qs_kind%all_potential)
    2147          480 :          CALL set_default_all_potential(qs_kind%all_potential, z, zeff_correction)
    2148          480 :          CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
    2149          480 :          IF (.NOT. ASSOCIATED(elec_conf)) THEN
    2150          480 :             CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
    2151          480 :             CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
    2152              :          END IF
    2153          480 :          CPASSERT(.NOT. qs_kind%floating)
    2154          480 :          IF (qs_kind%ghost) THEN
    2155            0 :             CALL get_qs_kind(qs_kind=qs_kind, elec_conf=elec_conf)
    2156            0 :             elec_conf(:) = 0
    2157              :             CALL get_potential(potential=qs_kind%all_potential, &
    2158            0 :                                elec_conf=elec_conf)
    2159            0 :             elec_conf(:) = 0
    2160              :             CALL set_potential(potential=qs_kind%all_potential, &
    2161              :                                zeff=0.0_dp, &
    2162            0 :                                zeff_correction=0.0_dp)
    2163              :          END IF
    2164              : 
    2165          480 :          check = ((potential_name /= '') .OR. explicit_potential) .AND. .NOT. silent
    2166              :          IF (check) &
    2167              :             CALL cp_warn(__LOCATION__, &
    2168              :                          "Information provided in the input file regarding POTENTIAL for KIND <"// &
    2169            0 :                          TRIM(qs_kind%name)//"> will be ignored!")
    2170              : 
    2171          480 :          check = ((k_rep > 0) .OR. explicit_basis) .AND. .NOT. silent
    2172              :          IF (check) &
    2173              :             CALL cp_warn(__LOCATION__, &
    2174              :                          "Information provided in the input file regarding BASIS for KIND <"// &
    2175           44 :                          TRIM(qs_kind%name)//"> will be ignored!")
    2176              : 
    2177              :       CASE (do_method_xtb)
    2178              :          ! Allocate all_potential
    2179         2208 :          CALL allocate_potential(qs_kind%all_potential)
    2180         2208 :          CALL set_default_all_potential(qs_kind%all_potential, z, zeff_correction)
    2181         2208 :          CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
    2182         2208 :          IF (.NOT. ASSOCIATED(elec_conf)) THEN
    2183         2208 :             CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
    2184         2208 :             CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
    2185              :          END IF
    2186         2208 :          CPASSERT(.NOT. qs_kind%floating)
    2187         2208 :          IF (qs_kind%ghost) THEN
    2188            0 :             CALL get_qs_kind(qs_kind=qs_kind, elec_conf=elec_conf)
    2189            0 :             elec_conf(:) = 0
    2190              :             CALL get_potential(potential=qs_kind%all_potential, &
    2191            0 :                                elec_conf=elec_conf)
    2192            0 :             elec_conf(:) = 0
    2193              :             CALL set_potential(potential=qs_kind%all_potential, &
    2194              :                                zeff=0.0_dp, &
    2195            0 :                                zeff_correction=0.0_dp)
    2196              :          END IF
    2197              : 
    2198         2208 :          check = ((potential_name /= '') .OR. explicit_potential) .AND. .NOT. silent
    2199              :          IF (check) &
    2200              :             CALL cp_warn(__LOCATION__, &
    2201              :                          "Information provided in the input file regarding POTENTIAL for KIND <"// &
    2202            0 :                          TRIM(qs_kind%name)//"> will be ignored!")
    2203              : 
    2204         2208 :          check = ((k_rep > 0) .OR. explicit_basis) .AND. .NOT. silent
    2205              :          IF (check) &
    2206              :             CALL cp_warn(__LOCATION__, &
    2207              :                          "Information provided in the input file regarding BASIS for KIND <"// &
    2208            0 :                          TRIM(qs_kind%name)//"> will be ignored!")
    2209              : 
    2210              :       CASE (do_method_pw)
    2211              :          ! PW DFT
    2212              :          ! Allocate and initialise the potential data set structure
    2213           22 :          IF (potential_name /= '') THEN
    2214           22 :             SELECT CASE (TRIM(potential_type))
    2215              :             CASE ("ALL", "ECP")
    2216              :                CALL cp_abort(__LOCATION__, &
    2217              :                              "PW DFT calculations only with potential type UPF or GTH possible."// &
    2218              :                              " <"//TRIM(potential_type)//"> was specified "// &
    2219            0 :                              "for the atomic kind <"//TRIM(qs_kind%name))
    2220              :             CASE ("GTH")
    2221            2 :                IF (potential_fn_kind == "-") THEN
    2222            2 :                   CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", c_val=potential_file_name)
    2223              :                ELSE
    2224            0 :                   potential_file_name = potential_fn_kind
    2225              :                END IF
    2226            2 :                CALL allocate_potential(qs_kind%gth_potential)
    2227              :                CALL read_potential(qs_kind%element_symbol, potential_name, &
    2228              :                                    qs_kind%gth_potential, zeff_correction, para_env, &
    2229            2 :                                    potential_file_name, potential_section, update_input)
    2230            2 :                CALL set_potential(qs_kind%gth_potential, z=z)
    2231            2 :                CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
    2232            2 :                IF (.NOT. ASSOCIATED(elec_conf)) THEN
    2233            2 :                   CALL get_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
    2234            2 :                   CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
    2235              :                ELSE
    2236            0 :                   CALL set_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
    2237              :                END IF
    2238              :             CASE ("UPF")
    2239         2100 :                ALLOCATE (qs_kind%upf_potential)
    2240           20 :                qs_kind%upf_potential%zion = 0
    2241           20 :                qs_kind%upf_potential%filename = ADJUSTL(TRIM(potential_name))
    2242           20 :                CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
    2243           20 :                IF (.NOT. ASSOCIATED(elec_conf)) THEN
    2244           20 :                   CALL set_qs_kind(qs_kind, elec_conf=qs_kind%upf_potential%econf)
    2245              :                END IF
    2246              :             CASE DEFAULT
    2247              :                CALL cp_abort(__LOCATION__, &
    2248              :                              "An invalid potential type <"// &
    2249              :                              TRIM(potential_type)//"> was specified "// &
    2250              :                              "for the atomic kind <"// &
    2251           22 :                              TRIM(qs_kind%name))
    2252              :             END SELECT
    2253              :          ELSE
    2254              :             CALL cp_abort(__LOCATION__, &
    2255              :                           "No potential type was defined for the "// &
    2256            0 :                           "atomic kind <"//TRIM(qs_kind%name)//">")
    2257              :          END IF
    2258              : 
    2259              :       CASE DEFAULT
    2260              : 
    2261              :          ! set ngauss for STO expansion
    2262         9447 :          CALL section_vals_val_get(dft_section, "QS%STO_NG", i_val=ngauss)
    2263              :          ! Allocate and initialise the basis set data set structure
    2264              :          ! first external basis sets
    2265        20805 :          DO i = 1, nb_rep
    2266        22712 :             SELECT CASE (basis_set_form(i))
    2267              :             CASE ("GTO")
    2268        11354 :                NULLIFY (tmp_basis_set)
    2269        11354 :                CALL allocate_gto_basis_set(tmp_basis_set)
    2270              :                CALL read_gto_basis_set(qs_kind%element_symbol, basis_set_name(i), &
    2271        11354 :                                        tmp_basis_set, para_env, dft_section)
    2272              :             CASE ("STO")
    2273            4 :                NULLIFY (sto_basis_set)
    2274            4 :                CALL allocate_sto_basis_set(sto_basis_set)
    2275              :                CALL read_sto_basis_set(qs_kind%element_symbol, basis_set_name(i), &
    2276            4 :                                        sto_basis_set, para_env, dft_section)
    2277            4 :                NULLIFY (tmp_basis_set)
    2278            4 :                CALL create_gto_from_sto_basis(sto_basis_set, tmp_basis_set, ngauss)
    2279            4 :                CALL deallocate_sto_basis_set(sto_basis_set)
    2280              :             CASE DEFAULT
    2281              :                CALL cp_abort(__LOCATION__, &
    2282              :                              "Invalid basis set form "//TRIM(basis_set_form(i))// &
    2283        11358 :                              "for atomic kind <"//TRIM(qs_kind%name)//">")
    2284              :             END SELECT
    2285        11358 :             tmp = basis_set_type(i)
    2286        11358 :             CALL uppercase(tmp)
    2287        20805 :             CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, tmp)
    2288              :          END DO
    2289              :          ! now explicit basis sets
    2290         9447 :          IF (explicit_basis) THEN
    2291          162 :             CALL section_vals_get(basis_section, n_repetition=nexp)
    2292          324 :             DO i = 1, nexp
    2293          162 :                NULLIFY (tmp_basis_set)
    2294          162 :                CALL allocate_gto_basis_set(tmp_basis_set)
    2295              :                CALL read_gto_basis_set(qs_kind%element_symbol, basis_type, &
    2296          162 :                                        tmp_basis_set, basis_section, i, dft_section)
    2297          162 :                tmp = basis_type
    2298          162 :                CALL uppercase(tmp)
    2299          324 :                CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, tmp)
    2300              :             END DO
    2301              :          END IF
    2302              :          ! combine multiple basis sets
    2303       198387 :          DO i = 1, SIZE(qs_kind%basis_sets)
    2304       188940 :             NULLIFY (tmp_basis_set)
    2305              :             CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
    2306       188940 :                                           inumbas=i, basis_type=basis_type)
    2307       188940 :             IF (basis_type == "") CYCLE
    2308        11520 :             jj = i
    2309       228197 :             DO j = i + 1, SIZE(qs_kind%basis_sets)
    2310       216677 :                jj = jj + 1
    2311       216677 :                NULLIFY (sup_basis_set)
    2312              :                CALL get_basis_from_container(qs_kind%basis_sets, basis_set=sup_basis_set, &
    2313       216677 :                                              inumbas=jj, basis_type=tmp)
    2314       228197 :                IF (basis_type == tmp) THEN
    2315              :                   ! we found a match, combine the basis sets and delete the second
    2316            0 :                   CALL combine_basis_sets(tmp_basis_set, sup_basis_set)
    2317            0 :                   CALL remove_basis_from_container(qs_kind%basis_sets, jj)
    2318            0 :                   jj = jj - 1
    2319              :                END IF
    2320              :             END DO
    2321       198387 :             NULLIFY (sup_basis_set)
    2322              :          END DO
    2323              : 
    2324              :          ! check that we have an orbital basis set
    2325         9447 :          nobasis = .TRUE.
    2326       198387 :          DO i = 1, SIZE(qs_kind%basis_sets)
    2327       188940 :             NULLIFY (tmp_basis_set)
    2328              :             CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
    2329       188940 :                                           inumbas=i, basis_type=basis_type)
    2330       198387 :             IF (basis_type == "ORB") nobasis = .FALSE.
    2331              :          END DO
    2332         9447 :          IF (nobasis) THEN
    2333              :             CALL cp_abort(__LOCATION__, &
    2334              :                           "No basis set type was defined for the "// &
    2335            0 :                           "atomic kind <"//TRIM(qs_kind%name)//">")
    2336              :          END IF
    2337              : 
    2338              :          ! If Ghost atom we don't need to allocate/initialize anything connected to POTENTIAL
    2339        26084 :          IF (qs_kind%ghost .OR. qs_kind%floating) THEN
    2340          150 :             IF (ASSOCIATED(qs_kind%elec_conf)) qs_kind%elec_conf = 0
    2341              :          ELSE
    2342              :             ! Allocate and initialise the potential data set structure
    2343         9297 :             IF ((potential_name /= '') .OR. explicit_potential) THEN
    2344              :                ! determine the pseudopotential file to search
    2345         9297 :                IF (potential_fn_kind == "-") THEN
    2346         9287 :                   CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", c_val=potential_file_name)
    2347              :                ELSE
    2348           10 :                   potential_file_name = potential_fn_kind
    2349              :                END IF
    2350              :                !
    2351        10353 :                SELECT CASE (TRIM(potential_type))
    2352              :                CASE ("ALL")
    2353         1056 :                   CALL allocate_potential(qs_kind%all_potential)
    2354              :                   CALL read_potential(qs_kind%element_symbol, potential_name, &
    2355              :                                       qs_kind%all_potential, zeff_correction, para_env, &
    2356         1056 :                                       potential_file_name, potential_section, update_input)
    2357         1056 :                   CALL set_potential(qs_kind%all_potential, z=z)
    2358         1056 :                   CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
    2359         1056 :                   IF (.NOT. ASSOCIATED(elec_conf)) THEN
    2360         1056 :                      CALL get_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
    2361         1056 :                      CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
    2362              :                   ELSE
    2363            0 :                      CALL set_potential(potential=qs_kind%all_potential, elec_conf=elec_conf)
    2364              :                   END IF
    2365              :                CASE ("GTH")
    2366         8201 :                   CALL allocate_potential(qs_kind%gth_potential)
    2367              :                   CALL read_potential(qs_kind%element_symbol, potential_name, &
    2368              :                                       qs_kind%gth_potential, zeff_correction, para_env, &
    2369         8201 :                                       potential_file_name, potential_section, update_input)
    2370         8201 :                   CALL set_potential(qs_kind%gth_potential, z=z)
    2371         8201 :                   CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
    2372         8201 :                   IF (.NOT. ASSOCIATED(elec_conf)) THEN
    2373         8197 :                      CALL get_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
    2374         8197 :                      CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
    2375              :                   ELSE
    2376            4 :                      CALL set_potential(potential=qs_kind%gth_potential, elec_conf=elec_conf)
    2377              :                   END IF
    2378              :                CASE ("ECP")
    2379           28 :                   CALL allocate_potential(qs_kind%sgp_potential)
    2380           28 :                   CALL get_potential(qs_kind%sgp_potential, description=description)
    2381              :                   CALL read_ecp_potential(ptable(z)%symbol, ecppot, &
    2382           28 :                                           potential_name, potential_file_name, potential_section)
    2383           28 :                   IF (ecp_semi_local) THEN
    2384           28 :                      description(1) = "Semi-local Gaussian pseudopotential                     "
    2385           28 :                      description(2) = "ECP "//TRIM(potential_name)
    2386           28 :                      description(3) = "LIBGRPP: A. V. Oleynichenko et al., Symmetry 15 197 2023"
    2387              :                      description(4) = "                                                        "
    2388              :                   ELSE
    2389            0 :                      description(4) = "ECP "//TRIM(potential_name)
    2390              :                   END IF
    2391              :                   CALL set_potential(qs_kind%sgp_potential, name=ecppot%pname, description=description, &
    2392              :                                      zeff=ecppot%zion, z=z, ecp_local=.TRUE., ecp_semi_local=ecp_semi_local, &
    2393              :                                      nloc=ecppot%nloc, nrloc=ecppot%nrloc, aloc=ecppot%aloc, bloc=ecppot%bloc, &
    2394           28 :                                      has_nlcc=.FALSE.)
    2395              :                   CALL set_potential(qs_kind%sgp_potential, sl_lmax=ecppot%lmax, &
    2396           28 :                                      npot=ecppot%npot, nrpot=ecppot%nrpot, apot=ecppot%apot, bpot=ecppot%bpot)
    2397              :                   ! convert PP
    2398           28 :                   IF (.NOT. ecp_semi_local) THEN
    2399            0 :                      CPABORT("ECPs are only well tested in their semi-local form")
    2400            0 :                      CALL get_qs_kind(qs_kind, basis_set=orb_basis_set)
    2401            0 :                      CALL sgp_construction(sgp_pot=sgppot, ecp_pot=ecppot, orb_basis=orb_basis_set, error=error)
    2402            0 :                      IF (iounit > 0 .AND. .NOT. silent) THEN
    2403            0 :                         WRITE (iounit, "(/,T2,'PP Transformation for ',A)") TRIM(ecppot%pname)
    2404            0 :                         IF (sgppot%has_local) THEN
    2405            0 :                            WRITE (iounit, "(T8,'Accuracy for local part:',T41,F10.3,'%',T61,F20.12)") error(4), error(1)
    2406              :                         END IF
    2407            0 :                         IF (sgppot%has_nonlocal) THEN
    2408            0 :                            WRITE (iounit, "(T8,'Accuracy for nonlocal part:',T41,F10.3,'%',T61,F20.12)") error(5), error(2)
    2409              :                         END IF
    2410            0 :                         IF (sgppot%has_nlcc) THEN
    2411            0 :                            WRITE (iounit, "(T8,'Accuracy for NLCC density:',T61,F20.12)") error(3)
    2412              :                         END IF
    2413              :                      END IF
    2414              :                   END IF
    2415           28 :                   IF (sgppot%has_nonlocal) THEN
    2416              :                      CALL set_potential(qs_kind%sgp_potential, n_nonlocal=sgppot%n_nonlocal, lmax=sgppot%lmax, &
    2417            0 :                                         is_nonlocal=sgppot%is_nonlocal)
    2418            0 :                      nnl = sgppot%n_nonlocal
    2419            0 :                      nppnl = 0
    2420            0 :                      DO l = 0, sgppot%lmax
    2421            0 :                         nppnl = nppnl + nnl*nco(l)
    2422              :                      END DO
    2423            0 :                      l = sgppot%lmax
    2424            0 :                      ALLOCATE (a_nl(nnl), h_nl(nnl, 0:l), c_nl(nnl, nnl, 0:l))
    2425            0 :                      a_nl(:) = sgppot%a_nonlocal(:)
    2426            0 :                      h_nl(:, :) = sgppot%h_nonlocal(:, :)
    2427            0 :                      DO l = 0, sgppot%lmax
    2428            0 :                         c_nl(:, :, l) = sgppot%c_nonlocal(:, :, l)*SQRT(2._dp*l + 1.0_dp)
    2429              :                      END DO
    2430            0 :                      CALL set_potential(qs_kind%sgp_potential, nppnl=nppnl, a_nonlocal=a_nl, h_nonlocal=h_nl, c_nonlocal=c_nl)
    2431              :                   ELSE
    2432           28 :                      CALL set_potential(qs_kind%sgp_potential, n_nonlocal=0, lmax=-1, is_nonlocal=sgppot%is_nonlocal)
    2433           28 :                      CALL set_potential(qs_kind%sgp_potential, nppnl=0)
    2434              :                   END IF
    2435              :                   !
    2436           28 :                   CPASSERT(.NOT. sgppot%has_local)
    2437           28 :                   CPASSERT(.NOT. sgppot%has_nlcc)
    2438              :                   ! core
    2439           28 :                   rc = 0.5_dp*qs_kind%covalent_radius*angstrom
    2440           28 :                   rc = MAX(rc, 0.2_dp)
    2441           28 :                   rc = MIN(rc, 1.0_dp)
    2442           28 :                   alpha = 1.0_dp/(2.0_dp*rc**2)
    2443           28 :                   ccore = ecppot%zion*SQRT((alpha/pi)**3)
    2444              :                   CALL set_potential(qs_kind%sgp_potential, alpha_core_charge=alpha, ccore_charge=ccore, &
    2445           28 :                                      core_charge_radius=rc)
    2446           28 :                   CALL atom_sgp_release(sgppot)
    2447           28 :                   CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
    2448           28 :                   IF (.NOT. ASSOCIATED(elec_conf)) THEN
    2449           28 :                      CALL set_qs_kind(qs_kind, elec_conf=ecppot%econf)
    2450              :                   END IF
    2451           28 :                   CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
    2452           28 :                   CALL set_potential(qs_kind%sgp_potential, elec_conf=elec_conf)
    2453              :                CASE ("UPF")
    2454           12 :                   CALL allocate_potential(qs_kind%sgp_potential)
    2455           12 :                   CALL get_potential(qs_kind%sgp_potential, description=description)
    2456           12 :                   description(4) = "UPF "//TRIM(potential_name)
    2457           12 :                   CALL atom_read_upf(upfpot, potential_name)
    2458              :                   CALL set_potential(qs_kind%sgp_potential, name=upfpot%pname, description=description, &
    2459           12 :                                      zeff=upfpot%zion, z=z, has_nlcc=upfpot%core_correction)
    2460              :                   ! convert pp
    2461           12 :                   CALL sgp_construction(sgp_pot=sgppot, upf_pot=upfpot, error=error)
    2462           12 :                   IF (iounit > 0 .AND. .NOT. silent) THEN
    2463            6 :                      WRITE (iounit, "(/,T2,'PP Transformation for ',A)") TRIM(upfpot%pname)
    2464            6 :                      IF (sgppot%has_local) THEN
    2465            6 :                         WRITE (iounit, "(T8,'Accuracy for local part:',T61,F20.12)") error(1)
    2466              :                      END IF
    2467            6 :                      IF (sgppot%has_nonlocal) THEN
    2468            3 :                         WRITE (iounit, "(T8,'Accuracy for nonlocal part:',T61,F20.12)") error(2)
    2469              :                      END IF
    2470            6 :                      IF (sgppot%has_nlcc) THEN
    2471            0 :                         WRITE (iounit, "(T8,'Accuracy for NLCC density:',T61,F20.12)") error(3)
    2472              :                      END IF
    2473              :                   END IF
    2474           12 :                   IF (sgppot%has_nonlocal) THEN
    2475              :                      CALL set_potential(qs_kind%sgp_potential, n_nonlocal=sgppot%n_nonlocal, lmax=sgppot%lmax, &
    2476            6 :                                         is_nonlocal=sgppot%is_nonlocal)
    2477            6 :                      nnl = sgppot%n_nonlocal
    2478            6 :                      nppnl = 0
    2479           12 :                      DO l = 0, sgppot%lmax
    2480           12 :                         nppnl = nppnl + nnl*nco(l)
    2481              :                      END DO
    2482            6 :                      l = sgppot%lmax
    2483           60 :                      ALLOCATE (a_nl(nnl), h_nl(nnl, 0:l), c_nl(nnl, nnl, 0:l))
    2484           54 :                      a_nl(:) = sgppot%a_nonlocal(:)
    2485           60 :                      h_nl(:, :) = sgppot%h_nonlocal(:, :)
    2486          444 :                      c_nl(:, :, :) = sgppot%c_nonlocal(:, :, :)
    2487            6 :                      CALL set_potential(qs_kind%sgp_potential, nppnl=nppnl, a_nonlocal=a_nl, h_nonlocal=h_nl, c_nonlocal=c_nl)
    2488              :                   ELSE
    2489            6 :                      CALL set_potential(qs_kind%sgp_potential, n_nonlocal=0, lmax=-1, is_nonlocal=sgppot%is_nonlocal)
    2490            6 :                      CALL set_potential(qs_kind%sgp_potential, nppnl=0)
    2491              :                   END IF
    2492           12 :                   CPASSERT(sgppot%has_local)
    2493              :                   ! core
    2494           12 :                   rc = sgppot%ac_local
    2495           12 :                   alpha = 1.0_dp/(2.0_dp*rc**2)
    2496           12 :                   ccore = upfpot%zion*SQRT((alpha/pi)**3)
    2497              :                   CALL set_potential(qs_kind%sgp_potential, alpha_core_charge=alpha, ccore_charge=ccore, &
    2498           12 :                                      core_charge_radius=rc)
    2499              :                   ! local potential
    2500           12 :                   nloc = sgppot%n_local
    2501           48 :                   ALLOCATE (aloc(nloc), cloc(nloc))
    2502          156 :                   aloc(1:nloc) = sgppot%a_local(1:nloc)
    2503          156 :                   cloc(1:nloc) = sgppot%c_local(1:nloc)
    2504           12 :                   CALL set_potential(qs_kind%sgp_potential, n_local=nloc, a_local=aloc, c_local=cloc)
    2505           12 :                   IF (sgppot%has_nlcc) THEN
    2506            0 :                      nlcc = sgppot%n_nlcc
    2507            0 :                      ALLOCATE (anlcc(nlcc), cnlcc(nlcc))
    2508            0 :                      anlcc(1:nlcc) = sgppot%a_nlcc(1:nlcc)
    2509            0 :                      cnlcc(1:nlcc) = sgppot%c_nlcc(1:nlcc)
    2510            0 :                      CALL set_potential(qs_kind%sgp_potential, has_nlcc=.TRUE., n_nlcc=nlcc, a_nlcc=anlcc, c_nlcc=cnlcc)
    2511              :                   END IF
    2512           12 :                   CALL set_potential(qs_kind%sgp_potential, z=z)
    2513           12 :                   CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
    2514           12 :                   IF (.NOT. ASSOCIATED(elec_conf)) THEN
    2515           12 :                      CALL set_qs_kind(qs_kind, elec_conf=upfpot%econf)
    2516              :                   END IF
    2517           12 :                   CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
    2518           12 :                   CALL set_potential(qs_kind%sgp_potential, elec_conf=elec_conf)
    2519           12 :                   CALL atom_release_upf(upfpot)
    2520           12 :                   CALL atom_sgp_release(sgppot)
    2521              :                CASE DEFAULT
    2522              :                   CALL cp_abort(__LOCATION__, &
    2523              :                                 "An invalid potential type <"// &
    2524              :                                 TRIM(potential_name)//"> was specified "// &
    2525              :                                 "for the atomic kind <"// &
    2526         9297 :                                 TRIM(qs_kind%name))
    2527              :                END SELECT
    2528              :             ELSE
    2529              :                CALL cp_abort(__LOCATION__, &
    2530              :                              "No potential type was defined for the "// &
    2531            0 :                              "atomic kind <"//TRIM(qs_kind%name)//">")
    2532              :             END IF
    2533              : 
    2534         9297 :             CALL check_potential_basis_compatibility(qs_kind)
    2535              : 
    2536              :             ! Allocate and initialise the potential data set structure
    2537         9297 :             IF ((kgpot_name /= '') .OR. explicit_kgpot) THEN
    2538         9297 :                ipos = INDEX(kgpot_name, "-")
    2539         9297 :                IF (ipos > 1) THEN
    2540           20 :                   kgpot_type = kgpot_name(:ipos - 1)
    2541              :                ELSE
    2542         9277 :                   kgpot_type = kgpot_name
    2543              :                END IF
    2544         9297 :                CALL uppercase(kgpot_type)
    2545              : 
    2546         9317 :                SELECT CASE (TRIM(kgpot_type))
    2547              :                CASE ("TNADD")
    2548              :                   ! determine the pseudopotential file to search
    2549           20 :                   IF (kg_potential_fn_kind == "-") THEN
    2550           20 :                      CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", c_val=potential_file_name)
    2551              :                   ELSE
    2552            0 :                      potential_file_name = kg_potential_fn_kind
    2553              :                   END IF
    2554           20 :                   CALL allocate_potential(qs_kind%tnadd_potential)
    2555              :                   CALL read_potential(qs_kind%element_symbol, kgpot_name, &
    2556              :                                       qs_kind%tnadd_potential, para_env, &
    2557           20 :                                       potential_file_name, kgpot_section, update_input)
    2558              :                CASE ("NONE")
    2559         9277 :                   NULLIFY (qs_kind%tnadd_potential)
    2560              :                CASE DEFAULT
    2561              :                   CALL cp_abort(__LOCATION__, &
    2562              :                                 "An invalid kg_potential type <"// &
    2563              :                                 TRIM(potential_name)//"> was specified "// &
    2564              :                                 "for the atomic kind <"// &
    2565         9297 :                                 TRIM(qs_kind%name))
    2566              :                END SELECT
    2567              :             END IF
    2568              :          END IF
    2569              :       END SELECT
    2570              : 
    2571        14397 :       CALL timestop(handle)
    2572              : 
    2573      8609406 :    END SUBROUTINE read_qs_kind
    2574              : 
    2575              : ! **************************************************************************************************
    2576              : !> \brief Ensure pseudo-potential and basis set were optimized for same number of valence electrons
    2577              : !> \param qs_kind ...
    2578              : !> \author Ole Schuett
    2579              : ! **************************************************************************************************
    2580         9297 :    SUBROUTINE check_potential_basis_compatibility(qs_kind)
    2581              :       TYPE(qs_kind_type), INTENT(INOUT)                  :: qs_kind
    2582              : 
    2583              :       CHARACTER(LEN=default_string_length)               :: name
    2584              :       INTEGER                                            :: nbs, npp
    2585              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
    2586              :       TYPE(gto_basis_set_type), POINTER                  :: basis_set
    2587              : 
    2588         9297 :       CALL get_qs_kind(qs_kind, name=name, gth_potential=gth_potential, basis_set=basis_set)
    2589              : 
    2590         9297 :       npp = -1; nbs = -1
    2591         9297 :       IF (ASSOCIATED(gth_potential)) &
    2592         8201 :          npp = parse_valence_electrons(gth_potential%aliases)
    2593         9297 :       IF (ASSOCIATED(basis_set)) &
    2594         9297 :          nbs = parse_valence_electrons(basis_set%aliases)
    2595              : 
    2596         9297 :       IF (npp >= 0 .AND. nbs >= 0 .AND. npp /= nbs) &
    2597              :          CALL cp_abort(__LOCATION__, "Basis-set and pseudo-potential of atomic kind '"//TRIM(name)//"'"// &
    2598            0 :                        " were optimized for different valence electron numbers.")
    2599              : 
    2600         9297 :    END SUBROUTINE check_potential_basis_compatibility
    2601              : 
    2602              : ! **************************************************************************************************
    2603              : !> \brief Tries to parse valence eletron number using "-QXXX" notation, returns -1 if not found.
    2604              : !> \param string ...
    2605              : !> \return ...
    2606              : !> \author Ole Schuett
    2607              : ! **************************************************************************************************
    2608        17498 :    FUNCTION parse_valence_electrons(string) RESULT(n)
    2609              :       CHARACTER(*)                                       :: string
    2610              :       INTEGER                                            :: n
    2611              : 
    2612              :       INTEGER                                            :: i, istat, j
    2613              : 
    2614        17498 :       i = INDEX(string, "-Q", .TRUE.)
    2615        17498 :       IF (i == 0) THEN
    2616         6078 :          n = -1
    2617              :       ELSE
    2618        11420 :          j = SCAN(string(i + 2:), "- ")
    2619        11420 :          READ (string(i + 2:i + j), '(I3)', iostat=istat) n
    2620        11420 :          IF (istat /= 0) n = -1
    2621              :       END IF
    2622              : 
    2623        17498 :    END FUNCTION
    2624              : 
    2625              : ! **************************************************************************************************
    2626              : !> \brief Read an atomic kind set data set from the input file.
    2627              : !> \param qs_kind_set ...
    2628              : !> \param atomic_kind_set ...
    2629              : !> \param kind_section ...
    2630              : !> \param para_env ...
    2631              : !> \param force_env_section ...
    2632              : !> \param silent ...
    2633              : ! **************************************************************************************************
    2634         7448 :    SUBROUTINE create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, para_env, &
    2635              :                                  force_env_section, silent)
    2636              : 
    2637              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2638              :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
    2639              :       TYPE(section_vals_type), POINTER                   :: kind_section
    2640              :       TYPE(mp_para_env_type), POINTER                    :: para_env
    2641              :       TYPE(section_vals_type), POINTER                   :: force_env_section
    2642              :       LOGICAL, INTENT(IN)                                :: silent
    2643              : 
    2644              :       CHARACTER(len=*), PARAMETER :: routineN = 'create_qs_kind_set'
    2645              : 
    2646              :       INTEGER                                            :: handle, ikind, method, nkind, qs_method
    2647              :       LOGICAL                                            :: no_fail
    2648              : 
    2649         7448 :       CALL timeset(routineN, handle)
    2650              : 
    2651         7448 :       IF (ASSOCIATED(qs_kind_set)) CPABORT("create_qs_kind_set: qs_kind_set already associated")
    2652         7448 :       IF (.NOT. ASSOCIATED(atomic_kind_set)) CPABORT("create_qs_kind_set: atomic_kind_set not associated")
    2653              : 
    2654         7448 :       no_fail = .FALSE.
    2655              : 
    2656              :       ! Between all methods only SE and DFTB/xTB may not need a KIND section.
    2657         7448 :       CALL section_vals_val_get(force_env_section, "METHOD", i_val=method)
    2658         7448 :       IF (method == do_qs) THEN
    2659         7428 :          CALL section_vals_val_get(force_env_section, "DFT%QS%METHOD", i_val=qs_method)
    2660          998 :          SELECT CASE (qs_method)
    2661              :          CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6fm, do_method_pm6, &
    2662              :                do_method_pdg, do_method_rm1, do_method_mndod, do_method_pnnl)
    2663          998 :             no_fail = .TRUE.
    2664              :          CASE (do_method_dftb)
    2665          222 :             no_fail = .TRUE.
    2666              :          CASE (do_method_xtb)
    2667         7428 :             no_fail = .TRUE.
    2668              :          END SELECT
    2669           20 :       ELSE IF (method == do_sirius) THEN
    2670           16 :          qs_method = do_method_pw
    2671              :       ELSE
    2672            4 :          qs_method = method
    2673              :       END IF
    2674              : 
    2675         7448 :       nkind = SIZE(atomic_kind_set)
    2676       185701 :       ALLOCATE (qs_kind_set(nkind))
    2677              : 
    2678        21845 :       DO ikind = 1, nkind
    2679        14397 :          qs_kind_set(ikind)%name = atomic_kind_set(ikind)%name
    2680        14397 :          qs_kind_set(ikind)%element_symbol = atomic_kind_set(ikind)%element_symbol
    2681        14397 :          qs_kind_set(ikind)%natom = atomic_kind_set(ikind)%natom
    2682              :          CALL read_qs_kind(qs_kind_set(ikind), kind_section, para_env, force_env_section, &
    2683        21845 :                            no_fail, qs_method, silent)
    2684              :       END DO
    2685              : 
    2686         7448 :       CALL timestop(handle)
    2687              : 
    2688        14896 :    END SUBROUTINE create_qs_kind_set
    2689              : 
    2690              : ! **************************************************************************************************
    2691              : !> \brief This routines should perform only checks. no settings are allowed at
    2692              : !>     this level anymore..
    2693              : !> \param qs_kind ...
    2694              : !> \param dft_control ...
    2695              : !> \param subsys_section ...
    2696              : ! **************************************************************************************************
    2697        14319 :    SUBROUTINE check_qs_kind(qs_kind, dft_control, subsys_section)
    2698              : 
    2699              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
    2700              :       TYPE(dft_control_type), INTENT(IN)                 :: dft_control
    2701              :       TYPE(section_vals_type), POINTER                   :: subsys_section
    2702              : 
    2703              :       INTEGER                                            :: gfn_type
    2704              :       LOGICAL                                            :: defined
    2705              :       TYPE(qs_dftb_atom_type), POINTER                   :: dftb_parameter
    2706              :       TYPE(semi_empirical_type), POINTER                 :: se_parameter
    2707              :       TYPE(xtb_atom_type), POINTER                       :: xtb_parameter
    2708              : 
    2709        14319 :       IF (dft_control%qs_control%semi_empirical) THEN
    2710         2240 :          CALL get_qs_kind(qs_kind, se_parameter=se_parameter)
    2711         2240 :          CPASSERT(ASSOCIATED(se_parameter))
    2712         2240 :          CALL get_se_param(se_parameter, defined=defined)
    2713         2240 :          CPASSERT(defined)
    2714         2240 :          CALL write_se_param(se_parameter, subsys_section)
    2715        12079 :       ELSE IF (dft_control%qs_control%dftb) THEN
    2716          480 :          CALL get_qs_kind(qs_kind, dftb_parameter=dftb_parameter)
    2717          480 :          CPASSERT(ASSOCIATED(dftb_parameter))
    2718          480 :          CALL get_dftb_atom_param(dftb_parameter, defined=defined)
    2719          480 :          CPASSERT(defined)
    2720          480 :          CALL write_dftb_atom_param(dftb_parameter, subsys_section)
    2721        11599 :       ELSE IF (dft_control%qs_control%xtb) THEN
    2722         2208 :          IF (.NOT. (dft_control%qs_control%xtb_control%do_tblite)) THEN
    2723         2096 :             CALL get_qs_kind(qs_kind, xtb_parameter=xtb_parameter)
    2724         2096 :             CPASSERT(ASSOCIATED(xtb_parameter))
    2725         2096 :             gfn_type = dft_control%qs_control%xtb_control%gfn_type
    2726         2096 :             CALL write_xtb_atom_param(xtb_parameter, gfn_type, subsys_section)
    2727              :          END IF
    2728              :       END IF
    2729              : 
    2730        14319 :    END SUBROUTINE check_qs_kind
    2731              : 
    2732              : ! **************************************************************************************************
    2733              : !> \brief ...
    2734              : !> \param qs_kind_set ...
    2735              : !> \param dft_control ...
    2736              : !> \param subsys_section ...
    2737              : ! **************************************************************************************************
    2738         7404 :    SUBROUTINE check_qs_kind_set(qs_kind_set, dft_control, subsys_section)
    2739              : 
    2740              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2741              :       TYPE(dft_control_type), INTENT(IN)                 :: dft_control
    2742              :       TYPE(section_vals_type), POINTER                   :: subsys_section
    2743              : 
    2744              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'check_qs_kind_set'
    2745              : 
    2746              :       INTEGER                                            :: handle, ikind, nkind
    2747              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
    2748              : 
    2749         7404 :       CALL timeset(routineN, handle)
    2750         7404 :       IF (ASSOCIATED(qs_kind_set)) THEN
    2751         7404 :          nkind = SIZE(qs_kind_set)
    2752        21723 :          DO ikind = 1, nkind
    2753        14319 :             qs_kind => qs_kind_set(ikind)
    2754        21723 :             CALL check_qs_kind(qs_kind, dft_control, subsys_section)
    2755              :          END DO
    2756         7404 :          IF (dft_control%qs_control%xtb) THEN
    2757              :             CALL write_xtb_kab_param(qs_kind_set, subsys_section, &
    2758          984 :                                      dft_control%qs_control%xtb_control)
    2759              :          END IF
    2760              :       ELSE
    2761            0 :          CPABORT("The pointer qs_kind_set is not associated")
    2762              :       END IF
    2763         7404 :       CALL timestop(handle)
    2764         7404 :    END SUBROUTINE check_qs_kind_set
    2765              : 
    2766              : ! **************************************************************************************************
    2767              : !> \brief ...
    2768              : !> \param qs_kind_set ...
    2769              : !> \param subsys_section ...
    2770              : !> \param xtb_control ...
    2771              : ! **************************************************************************************************
    2772          984 :    SUBROUTINE write_xtb_kab_param(qs_kind_set, subsys_section, xtb_control)
    2773              : 
    2774              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    2775              :       TYPE(section_vals_type), POINTER                   :: subsys_section
    2776              :       TYPE(xtb_control_type), POINTER                    :: xtb_control
    2777              : 
    2778              :       CHARACTER(LEN=default_string_length)               :: aname, bname
    2779              :       INTEGER                                            :: ikind, io_unit, jkind, nkind, za, zb
    2780              :       TYPE(cp_logger_type), POINTER                      :: logger
    2781              :       TYPE(qs_kind_type), POINTER                        :: qs_kinda, qs_kindb
    2782              :       TYPE(xtb_atom_type), POINTER                       :: xtb_parameter_a, xtb_parameter_b
    2783              : 
    2784          984 :       NULLIFY (logger)
    2785          984 :       logger => cp_get_default_logger()
    2786          984 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, subsys_section, &
    2787              :                                            "PRINT%KINDS/POTENTIAL"), cp_p_file)) THEN
    2788              : 
    2789            0 :          io_unit = cp_print_key_unit_nr(logger, subsys_section, "PRINT%KINDS", extension=".Log")
    2790            0 :          IF (io_unit > 0) THEN
    2791              : 
    2792            0 :             WRITE (io_unit, "(/,T2,A)") "xTB| Kab parameters"
    2793            0 :             nkind = SIZE(qs_kind_set)
    2794            0 :             DO ikind = 1, nkind
    2795            0 :                qs_kinda => qs_kind_set(ikind)
    2796            0 :                CALL get_qs_kind(qs_kinda, xtb_parameter=xtb_parameter_a)
    2797            0 :                CALL get_xtb_atom_param(xtb_parameter_a, aname=aname, z=za)
    2798            0 :                DO jkind = ikind, nkind
    2799            0 :                   qs_kindb => qs_kind_set(jkind)
    2800            0 :                   CALL get_qs_kind(qs_kindb, xtb_parameter=xtb_parameter_b)
    2801            0 :                   CALL get_xtb_atom_param(xtb_parameter_b, aname=bname, z=zb)
    2802              :                   WRITE (io_unit, "(A,T10,A15,T25,A15,T71,F10.3)") &
    2803            0 :                      "    Kab:", TRIM(aname), TRIM(bname), xtb_set_kab(za, zb, xtb_control)
    2804              :                END DO
    2805              :             END DO
    2806            0 :             WRITE (io_unit, *)
    2807              : 
    2808              :          END IF
    2809              : 
    2810            0 :          CALL cp_print_key_finished_output(io_unit, logger, subsys_section, "PRINT%KINDS")
    2811              :       END IF
    2812              : 
    2813          984 :    END SUBROUTINE write_xtb_kab_param
    2814              : 
    2815              : ! **************************************************************************************************
    2816              : !> \brief Set the components of an atomic kind data set.
    2817              : !> \param qs_kind ...
    2818              : !> \param paw_atom ...
    2819              : !> \param ghost ...
    2820              : !> \param floating ...
    2821              : !> \param hard_radius ...
    2822              : !> \param hard0_radius ...
    2823              : !> \param covalent_radius ...
    2824              : !> \param vdw_radius ...
    2825              : !> \param lmax_rho0 ...
    2826              : !> \param zeff ...
    2827              : !> \param no_optimize ...
    2828              : !> \param dispersion ...
    2829              : !> \param u_minus_j ...
    2830              : !> \param reltmat ...
    2831              : !> \param dftb_parameter ...
    2832              : !> \param xtb_parameter ...
    2833              : !> \param elec_conf ...
    2834              : !> \param pao_basis_size ...
    2835              : ! **************************************************************************************************
    2836        23417 :    SUBROUTINE set_qs_kind(qs_kind, paw_atom, ghost, floating, hard_radius, hard0_radius, &
    2837              :                           covalent_radius, vdw_radius, lmax_rho0, zeff, &
    2838              :                           no_optimize, dispersion, u_minus_j, reltmat, &
    2839              :                           dftb_parameter, xtb_parameter, &
    2840        23417 :                           elec_conf, pao_basis_size)
    2841              : 
    2842              :       TYPE(qs_kind_type), INTENT(INOUT)                  :: qs_kind
    2843              :       LOGICAL, INTENT(IN), OPTIONAL                      :: paw_atom, ghost, floating
    2844              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: hard_radius, hard0_radius, &
    2845              :                                                             covalent_radius, vdw_radius
    2846              :       INTEGER, INTENT(IN), OPTIONAL                      :: lmax_rho0
    2847              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: zeff
    2848              :       LOGICAL, INTENT(IN), OPTIONAL                      :: no_optimize
    2849              :       TYPE(qs_atom_dispersion_type), OPTIONAL, POINTER   :: dispersion
    2850              :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: u_minus_j
    2851              :       REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER  :: reltmat
    2852              :       TYPE(qs_dftb_atom_type), OPTIONAL, POINTER         :: dftb_parameter
    2853              :       TYPE(xtb_atom_type), OPTIONAL, POINTER             :: xtb_parameter
    2854              :       INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL        :: elec_conf
    2855              :       INTEGER, INTENT(IN), OPTIONAL                      :: pao_basis_size
    2856              : 
    2857        23417 :       IF (PRESENT(dftb_parameter)) qs_kind%dftb_parameter => dftb_parameter
    2858        23417 :       IF (PRESENT(xtb_parameter)) qs_kind%xtb_parameter => xtb_parameter
    2859        23417 :       IF (PRESENT(elec_conf)) THEN
    2860        14247 :          IF (ASSOCIATED(qs_kind%elec_conf)) THEN
    2861            0 :             DEALLOCATE (qs_kind%elec_conf)
    2862              :          END IF
    2863        42741 :          ALLOCATE (qs_kind%elec_conf(0:SIZE(elec_conf) - 1))
    2864        51799 :          qs_kind%elec_conf(:) = elec_conf(:)
    2865              :       END IF
    2866        23417 :       IF (PRESENT(paw_atom)) qs_kind%paw_atom = paw_atom
    2867        23417 :       IF (PRESENT(hard_radius)) qs_kind%hard_radius = hard_radius
    2868        23417 :       IF (PRESENT(hard0_radius)) qs_kind%hard0_radius = hard0_radius
    2869        23417 :       IF (PRESENT(covalent_radius)) qs_kind%covalent_radius = covalent_radius
    2870        23417 :       IF (PRESENT(vdw_radius)) qs_kind%vdw_radius = vdw_radius
    2871        23417 :       IF (PRESENT(lmax_rho0)) qs_kind%lmax_rho0 = lmax_rho0
    2872        23417 :       IF (PRESENT(zeff)) THEN
    2873            0 :          IF (ASSOCIATED(qs_kind%all_potential)) THEN
    2874            0 :             CALL set_potential(potential=qs_kind%all_potential, zeff=zeff)
    2875            0 :          ELSE IF (ASSOCIATED(qs_kind%gth_potential)) THEN
    2876            0 :             CALL set_potential(potential=qs_kind%gth_potential, zeff=zeff)
    2877            0 :          ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
    2878            0 :             CALL set_potential(potential=qs_kind%sgp_potential, zeff=zeff)
    2879              :          END IF
    2880              :       END IF
    2881        23417 :       IF (PRESENT(ghost)) qs_kind%ghost = ghost
    2882              : 
    2883        23417 :       IF (PRESENT(floating)) qs_kind%floating = floating
    2884              : 
    2885        23417 :       IF (PRESENT(no_optimize)) qs_kind%no_optimize = no_optimize
    2886              : 
    2887        23417 :       IF (PRESENT(dispersion)) qs_kind%dispersion => dispersion
    2888              : 
    2889        23417 :       IF (PRESENT(u_minus_j)) THEN
    2890          476 :          IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
    2891          476 :             qs_kind%dft_plus_u%u_minus_j = u_minus_j
    2892              :          END IF
    2893              :       END IF
    2894              : 
    2895        23417 :       IF (PRESENT(reltmat)) qs_kind%reltmat => reltmat
    2896              : 
    2897        23417 :       IF (PRESENT(pao_basis_size)) qs_kind%pao_basis_size = pao_basis_size
    2898              : 
    2899        23417 :    END SUBROUTINE set_qs_kind
    2900              : 
    2901              : ! **************************************************************************************************
    2902              : !> \brief Write an atomic kind data set to the output unit.
    2903              : !> \param qs_kind ...
    2904              : !> \param kind_number ...
    2905              : !> \param output_unit ...
    2906              : !> \par History
    2907              : !>      Creation (09.02.2002,MK)
    2908              : ! **************************************************************************************************
    2909         3632 :    SUBROUTINE write_qs_kind(qs_kind, kind_number, output_unit)
    2910              : 
    2911              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
    2912              :       INTEGER, INTENT(in)                                :: kind_number, output_unit
    2913              : 
    2914              :       CHARACTER(LEN=3)                                   :: yon
    2915              :       CHARACTER(LEN=default_string_length)               :: basis_type, bstring
    2916              :       INTEGER                                            :: ibas
    2917              :       LOGICAL                                            :: do_print
    2918              :       TYPE(gto_basis_set_type), POINTER                  :: tmp_basis
    2919              : 
    2920         3632 :       IF (output_unit > 0) THEN
    2921              : 
    2922         3632 :          IF (ASSOCIATED(qs_kind)) THEN
    2923              :             WRITE (UNIT=output_unit, FMT="(/,T2,I2,A,T57,A,T75,I6)") &
    2924         3632 :                kind_number, ". Atomic kind: "//TRIM(qs_kind%name), &
    2925         7264 :                "Number of atoms: ", qs_kind%natom
    2926              : 
    2927        76272 :             DO ibas = 1, SIZE(qs_kind%basis_sets, 1)
    2928        72640 :                NULLIFY (tmp_basis)
    2929              :                CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
    2930        72640 :                                              inumbas=ibas, basis_type=basis_type)
    2931        72640 :                do_print = .TRUE.
    2932        67965 :                SELECT CASE (basis_type)
    2933              :                CASE DEFAULT
    2934        67965 :                   bstring = "Basis Set"
    2935         3546 :                   do_print = .FALSE.
    2936              :                CASE ("ORB")
    2937         3546 :                   bstring = "Orbital Basis Set"
    2938              :                CASE ("ORB_SOFT")
    2939          466 :                   bstring = "GAPW Soft Basis Set"
    2940            0 :                   do_print = .FALSE.
    2941              :                CASE ("AUX")
    2942            0 :                   bstring = "Auxiliary Basis Set"
    2943              :                CASE ("MIN")
    2944            0 :                   bstring = "Minimal Basis Set"
    2945              :                CASE ("RI_AUX")
    2946          344 :                   bstring = "RI Auxiliary Basis Set"
    2947              :                CASE ("AUX_FIT")
    2948          220 :                   bstring = "Auxiliary Fit Basis Set"
    2949              :                CASE ("LRI_AUX")
    2950           15 :                   bstring = "LRI Basis Set"
    2951              :                CASE ("P_LRI_AUX")
    2952            4 :                   bstring = "LRI Basis Set for TDDFPT"
    2953              :                CASE ("RI_XAS")
    2954            0 :                   bstring = "RI XAS Basis Set"
    2955              :                CASE ("RI_HFX")
    2956        72640 :                   bstring = "RI HFX Basis Set"
    2957              :                END SELECT
    2958              : 
    2959         3632 :                IF (do_print) THEN
    2960         4209 :                   CALL write_orb_basis_set(tmp_basis, output_unit, bstring)
    2961              :                END IF
    2962              : 
    2963              :             END DO
    2964              : 
    2965         3632 :             IF (qs_kind%ghost) THEN
    2966              :                WRITE (UNIT=output_unit, FMT="(/,T6,A)") &
    2967           11 :                   "The atoms of this atomic kind are GHOST atoms!"
    2968              :             END IF
    2969         3632 :             IF (qs_kind%floating) THEN
    2970              :                WRITE (UNIT=output_unit, FMT="(/,T6,A)") &
    2971            0 :                   "The atoms of this atomic kind are FLOATING BASIS FUNCTIONS."
    2972              :             END IF
    2973         3632 :             IF (qs_kind%covalent_radius > 0.0_dp) THEN
    2974              :                WRITE (UNIT=output_unit, FMT="(/,T8,A,T71,F10.3)") &
    2975         2426 :                   "Atomic covalent radius [Angstrom]:", &
    2976         4852 :                   qs_kind%covalent_radius*angstrom
    2977              :             END IF
    2978         3632 :             IF (qs_kind%vdw_radius > 0.0_dp) THEN
    2979              :                WRITE (UNIT=output_unit, FMT="(/,T8,A,T71,F10.3)") &
    2980         2426 :                   "Atomic van der Waals radius [Angstrom]:", &
    2981         4852 :                   qs_kind%vdw_radius*angstrom
    2982              :             END IF
    2983         3632 :             IF (qs_kind%paw_atom) THEN
    2984              :                WRITE (UNIT=output_unit, FMT="(/,T6,A)") &
    2985          376 :                   "The atoms of this atomic kind are PAW atoms (GAPW):"
    2986              :                WRITE (UNIT=output_unit, FMT="(T8,A,T71,F10.3)") &
    2987          376 :                   "Hard Gaussian function radius:", qs_kind%hard_radius, &
    2988          376 :                   "Rho0 radius:", qs_kind%hard0_radius, &
    2989          376 :                   "Maximum GTO radius used for PAW projector construction:", &
    2990          752 :                   qs_kind%max_rad_local
    2991          376 :                NULLIFY (tmp_basis)
    2992              :                CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
    2993          376 :                                              basis_type="ORB_SOFT")
    2994          376 :                CALL write_orb_basis_set(tmp_basis, output_unit, "GAPW Soft Basis Set")
    2995              :             END IF
    2996              :             ! Potentials
    2997         3632 :             IF (ASSOCIATED(qs_kind%all_potential)) CALL write_potential(qs_kind%all_potential, output_unit)
    2998         3632 :             IF (ASSOCIATED(qs_kind%gth_potential)) CALL write_potential(qs_kind%gth_potential, output_unit)
    2999         3632 :             IF (ASSOCIATED(qs_kind%sgp_potential)) CALL write_potential(qs_kind%sgp_potential, output_unit)
    3000         3632 :             IF (ASSOCIATED(qs_kind%tnadd_potential)) CALL write_potential(qs_kind%tnadd_potential, output_unit)
    3001         3632 :             IF (ASSOCIATED(qs_kind%dft_plus_u)) THEN
    3002              :                WRITE (UNIT=output_unit, FMT="(/,T6,A,/,T8,A,T76,I5,/,T8,A,T73,F8.3)") &
    3003           16 :                   "A DFT+U correction is applied to atoms of this atomic kind:", &
    3004           16 :                   "Angular quantum momentum number L:", qs_kind%dft_plus_u%l, &
    3005           32 :                   "U(eff) = (U - J) value in [eV]:", qs_kind%dft_plus_u%u_minus_j_target*evolt
    3006           16 :                IF (qs_kind%dft_plus_u%u_ramping > 0.0_dp) THEN
    3007            4 :                   IF (qs_kind%dft_plus_u%init_u_ramping_each_scf) THEN
    3008            2 :                      yon = "YES"
    3009              :                   ELSE
    3010            2 :                      yon = " NO"
    3011              :                   END IF
    3012              :                   WRITE (UNIT=output_unit, FMT="(T8,A,T73,F8.3,/,T8,A,T73,ES8.1,/,T8,A,T78,A3)") &
    3013            4 :                      "Increment for U ramping in [eV]:", qs_kind%dft_plus_u%u_ramping*evolt, &
    3014            4 :                      "SCF threshold value for U ramping:", qs_kind%dft_plus_u%eps_u_ramping, &
    3015            8 :                      "Set U ramping value to zero before each wavefunction optimisation:", yon
    3016              :                END IF
    3017           16 :                IF (ASSOCIATED(qs_kind%dft_plus_u%orbitals)) THEN
    3018              :                   WRITE (UNIT=output_unit, FMT="(T8,A)") &
    3019            2 :                      "An initial orbital occupation is requested:"
    3020            2 :                   IF (ASSOCIATED(qs_kind%dft_plus_u%nelec)) THEN
    3021            4 :                      IF (ANY(qs_kind%dft_plus_u%nelec(:) >= 0.5_dp)) THEN
    3022            0 :                         IF (SIZE(qs_kind%dft_plus_u%nelec) > 1) THEN
    3023              :                            WRITE (UNIT=output_unit, FMT="(T9,A,T75,F6.2)") &
    3024            0 :                               "Number of alpha electrons:", &
    3025            0 :                               qs_kind%dft_plus_u%nelec(1), &
    3026            0 :                               "Number of beta electrons:", &
    3027            0 :                               qs_kind%dft_plus_u%nelec(2)
    3028              :                         ELSE
    3029              :                            WRITE (UNIT=output_unit, FMT="(T9,A,T75,F6.2)") &
    3030            0 :                               "Number of electrons:", &
    3031            0 :                               qs_kind%dft_plus_u%nelec(1)
    3032              :                         END IF
    3033              :                      END IF
    3034              :                   END IF
    3035              :                   WRITE (UNIT=output_unit, FMT="(T9,A,(T78,I3))") &
    3036            2 :                      "Preferred (initial) orbital occupation order (orbital M values):", &
    3037            4 :                      qs_kind%dft_plus_u%orbitals(:)
    3038              :                   WRITE (UNIT=output_unit, FMT="(T9,A,T71,ES10.3,/,T9,A,T76,I5)") &
    3039            2 :                      "Threshold value for the SCF convergence criterion:", &
    3040            2 :                      qs_kind%dft_plus_u%eps_scf, &
    3041            2 :                      "Number of initial SCF iterations:", &
    3042            4 :                      qs_kind%dft_plus_u%max_scf
    3043            2 :                   IF (qs_kind%dft_plus_u%smear) THEN
    3044              :                      WRITE (UNIT=output_unit, FMT="(T9,A)") &
    3045            2 :                         "A smearing of the orbital occupations will be performed"
    3046              :                   END IF
    3047              :                END IF
    3048              :             END IF
    3049              :          ELSE
    3050            0 :             CPABORT("")
    3051              :          END IF
    3052              : 
    3053              :       END IF
    3054              : 
    3055         3632 :    END SUBROUTINE write_qs_kind
    3056              : 
    3057              : ! **************************************************************************************************
    3058              : !> \brief Write an atomic kind set data set to the output unit.
    3059              : !> \param qs_kind_set ...
    3060              : !> \param subsys_section ...
    3061              : !> \par History
    3062              : !>      Creation (09.02.2002,MK)
    3063              : ! **************************************************************************************************
    3064         7414 :    SUBROUTINE write_qs_kind_set(qs_kind_set, subsys_section)
    3065              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    3066              :       TYPE(section_vals_type), POINTER                   :: subsys_section
    3067              : 
    3068              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'write_qs_kind_set'
    3069              : 
    3070              :       INTEGER                                            :: handle, ikind, nkind, output_unit
    3071              :       TYPE(cp_logger_type), POINTER                      :: logger
    3072              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
    3073              : 
    3074         7414 :       CALL timeset(routineN, handle)
    3075              : 
    3076         7414 :       NULLIFY (logger)
    3077         7414 :       logger => cp_get_default_logger()
    3078              :       output_unit = cp_print_key_unit_nr(logger, subsys_section, &
    3079         7414 :                                          "PRINT%KINDS", extension=".Log")
    3080         7414 :       IF (output_unit > 0) THEN
    3081         1930 :          IF (ASSOCIATED(qs_kind_set)) THEN
    3082         1930 :             WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") "ATOMIC KIND INFORMATION"
    3083         1930 :             nkind = SIZE(qs_kind_set)
    3084         5562 :             DO ikind = 1, nkind
    3085         3632 :                qs_kind => qs_kind_set(ikind)
    3086         5562 :                CALL write_qs_kind(qs_kind, ikind, output_unit)
    3087              :             END DO
    3088              :          ELSE
    3089            0 :             CPABORT("")
    3090              :          END IF
    3091              :       END IF
    3092              : 
    3093              :       CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
    3094         7414 :                                         "PRINT%KINDS")
    3095              : 
    3096         7414 :       CALL timestop(handle)
    3097              : 
    3098         7414 :    END SUBROUTINE write_qs_kind_set
    3099              : 
    3100              : ! **************************************************************************************************
    3101              : !> \brief Write all the GTO basis sets of an atomic kind set to the output
    3102              : !>     unit (for the printing of the unnormalized basis sets as read from
    3103              : !>           database).
    3104              : !> \param qs_kind_set ...
    3105              : !> \param subsys_section ...
    3106              : !> \par History
    3107              : !>      Creation (17.01.2002,MK)
    3108              : ! **************************************************************************************************
    3109         7398 :    SUBROUTINE write_gto_basis_sets(qs_kind_set, subsys_section)
    3110              : 
    3111              :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
    3112              :       TYPE(section_vals_type), POINTER                   :: subsys_section
    3113              : 
    3114              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'write_gto_basis_sets'
    3115              : 
    3116              :       CHARACTER(LEN=default_string_length)               :: basis_type, bstring
    3117              :       INTEGER                                            :: handle, ibas, ikind, nkind, output_unit
    3118              :       TYPE(cp_logger_type), POINTER                      :: logger
    3119              :       TYPE(gto_basis_set_type), POINTER                  :: tmp_basis
    3120              :       TYPE(qs_kind_type), POINTER                        :: qs_kind
    3121              : 
    3122         7398 :       CALL timeset(routineN, handle)
    3123              : 
    3124         7398 :       NULLIFY (logger)
    3125         7398 :       logger => cp_get_default_logger()
    3126              :       output_unit = cp_print_key_unit_nr(logger, subsys_section, &
    3127              :                                          "PRINT%KINDS/BASIS_SET", &
    3128         7398 :                                          extension=".Log")
    3129         7398 :       IF (output_unit > 0) THEN
    3130           60 :          IF (ASSOCIATED(qs_kind_set)) THEN
    3131              :             WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
    3132           60 :                "BASIS SET INFORMATION (Unnormalised Gaussian-type functions)"
    3133           60 :             nkind = SIZE(qs_kind_set)
    3134          175 :             DO ikind = 1, nkind
    3135          115 :                qs_kind => qs_kind_set(ikind)
    3136              :                WRITE (UNIT=output_unit, FMT="(/,T2,I2,A)") &
    3137          115 :                   ikind, ". Atomic kind: "//TRIM(qs_kind%name)
    3138              : 
    3139         2475 :                DO ibas = 1, SIZE(qs_kind%basis_sets, 1)
    3140         2300 :                   NULLIFY (tmp_basis)
    3141              :                   CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
    3142         2300 :                                                 inumbas=ibas, basis_type=basis_type)
    3143         2300 :                   IF (basis_type == "") CYCLE
    3144           11 :                   SELECT CASE (basis_type)
    3145              :                   CASE DEFAULT
    3146           11 :                      bstring = "Basis Set"
    3147              :                   CASE ("ORB")
    3148          115 :                      bstring = "Orbital Basis Set"
    3149              :                   CASE ("ORB_SOFT")
    3150           11 :                      bstring = "GAPW Soft Basis Set"
    3151              :                   CASE ("AUX")
    3152            0 :                      bstring = "Auxiliary Basis Set"
    3153              :                   CASE ("MIN")
    3154            0 :                      bstring = "Minimal Basis Set"
    3155              :                   CASE ("RI_AUX")
    3156            0 :                      bstring = "RI Auxiliary Basis Set"
    3157              :                   CASE ("AUX_FIT")
    3158            0 :                      bstring = "Auxiliary Fit Basis Set"
    3159              :                   CASE ("LRI_AUX")
    3160            2 :                      bstring = "LRI Basis Set"
    3161              :                   CASE ("P_LRI_AUX")
    3162            0 :                      bstring = "LRI Basis Set for TDDFPT"
    3163              :                   CASE ("RI_HFX")
    3164          139 :                      bstring = "RI HFX Basis Set"
    3165              :                   END SELECT
    3166              : 
    3167          254 :                   IF (ASSOCIATED(tmp_basis)) CALL write_gto_basis_set(tmp_basis, output_unit, bstring)
    3168              : 
    3169              :                END DO
    3170              : 
    3171              :             END DO
    3172              :          ELSE
    3173            0 :             CPABORT("")
    3174              :          END IF
    3175              :       END IF
    3176              : 
    3177              :       CALL cp_print_key_finished_output(output_unit, logger, subsys_section, &
    3178         7398 :                                         "PRINT%KINDS/BASIS_SET")
    3179              : 
    3180         7398 :       CALL timestop(handle)
    3181              : 
    3182         7398 :    END SUBROUTINE write_gto_basis_sets
    3183              : 
    3184              : ! **************************************************************************************************
    3185              : !> \brief ...
    3186              : !> \param atomic_kind ...
    3187              : !> \param qs_kind ...
    3188              : !> \param ncalc ...
    3189              : !> \param ncore ...
    3190              : !> \param nelem ...
    3191              : !> \param edelta ...
    3192              : ! **************************************************************************************************
    3193        90594 :    SUBROUTINE init_atom_electronic_state(atomic_kind, qs_kind, ncalc, ncore, nelem, edelta)
    3194              : 
    3195              :       TYPE(atomic_kind_type), INTENT(IN)                 :: atomic_kind
    3196              :       TYPE(qs_kind_type), INTENT(IN)                     :: qs_kind
    3197              :       INTEGER, DIMENSION(0:lmat, 10), INTENT(OUT)        :: ncalc, ncore, nelem
    3198              :       REAL(KIND=dp), DIMENSION(0:lmat, 10, 2), &
    3199              :          INTENT(OUT)                                     :: edelta
    3200              : 
    3201              :       INTEGER                                            :: i, ii, is, l, ll, ne, nn, z
    3202        45297 :       INTEGER, DIMENSION(:), POINTER                     :: econf
    3203        45297 :       INTEGER, DIMENSION(:, :), POINTER                  :: addel, laddel, naddel
    3204              :       LOGICAL                                            :: bs_occupation
    3205              :       REAL(KIND=dp)                                      :: dmag, magnetization
    3206              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
    3207              :       TYPE(sgp_potential_type), POINTER                  :: sgp_potential
    3208              : 
    3209        45297 :       CALL get_atomic_kind(atomic_kind, z=z)
    3210        45297 :       NULLIFY (gth_potential)
    3211              :       CALL get_qs_kind(qs_kind, &
    3212              :                        gth_potential=gth_potential, &
    3213              :                        sgp_potential=sgp_potential, &
    3214              :                        magnetization=magnetization, &
    3215              :                        bs_occupation=bs_occupation, &
    3216        45297 :                        addel=addel, laddel=laddel, naddel=naddel)
    3217              : 
    3218              :       ! electronic state
    3219        45297 :       nelem = 0
    3220        45297 :       ncore = 0
    3221        45297 :       ncalc = 0
    3222        45297 :       edelta = 0.0_dp
    3223        45297 :       IF (ASSOCIATED(gth_potential)) THEN
    3224        24165 :          CALL get_potential(gth_potential, elec_conf=econf)
    3225        24165 :          CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
    3226        21132 :       ELSE IF (ASSOCIATED(sgp_potential)) THEN
    3227           90 :          CALL get_potential(sgp_potential, elec_conf=econf)
    3228           90 :          CALL set_pseudo_state(econf, z, ncalc, ncore, nelem)
    3229              :       ELSE
    3230       105210 :          DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
    3231        84168 :             ll = 2*(2*l + 1)
    3232        84168 :             nn = ptable(z)%e_conv(l)
    3233        84168 :             ii = 0
    3234        21042 :             DO
    3235       120234 :                ii = ii + 1
    3236       120234 :                IF (nn <= ll) THEN
    3237        84168 :                   nelem(l, ii) = nn
    3238              :                   EXIT
    3239              :                ELSE
    3240        36066 :                   nelem(l, ii) = ll
    3241        36066 :                   nn = nn - ll
    3242              :                END IF
    3243              :             END DO
    3244              :          END DO
    3245      1493982 :          ncalc = nelem - ncore
    3246              :       END IF
    3247              : 
    3248              :       ! readjust the occupation number of the orbitals as requested by user
    3249              :       ! this is done to break symmetry (bs) and bias the initial guess
    3250              :       ! to the pre-defined multiplicity/charge state of the atom
    3251        45297 :       IF (bs_occupation) THEN
    3252          648 :          DO is = 1, 2
    3253         1176 :             DO i = 1, SIZE(addel, 1)
    3254          528 :                ne = addel(i, is)
    3255          528 :                l = laddel(i, is)
    3256          528 :                nn = naddel(i, is) - l
    3257          960 :                IF (ne /= 0) THEN
    3258          500 :                   IF (nn == 0) THEN
    3259            0 :                      DO ii = SIZE(nelem, 2), 1, -1
    3260            0 :                         IF (ncalc(l, ii) > 0) THEN
    3261            0 :                            IF ((ncalc(l, ii) + ne) < 2*(2*l + 1) + 1) THEN
    3262            0 :                               edelta(l, ii, is) = edelta(l, ii, is) + ne
    3263            0 :                               nn = ii
    3264              :                            ELSE
    3265            0 :                               edelta(l, ii + 1, is) = edelta(l, ii + 1, is) + ne
    3266            0 :                               nn = ii + 1
    3267              :                            END IF
    3268              :                            EXIT
    3269            0 :                         ELSE IF (ii == 1) THEN
    3270            0 :                            edelta(l, ii, is) = edelta(l, ii, is) + ne
    3271            0 :                            nn = ii
    3272              :                         END IF
    3273              :                      END DO
    3274              :                   ELSE
    3275          500 :                      edelta(l, nn, is) = edelta(l, nn, is) + ne
    3276              :                   END IF
    3277          500 :                   IF (ncalc(l, nn) + edelta(l, nn, is) < 0) THEN
    3278            0 :                      edelta(l, nn, is) = -ncalc(l, nn)
    3279              :                   END IF
    3280              :                END IF
    3281              :             END DO
    3282              :          END DO
    3283        30888 :          edelta = 0.5_dp*edelta
    3284        45081 :       ELSE IF (magnetization /= 0.0_dp) THEN
    3285            0 :          dmag = 0.5_dp*ABS(magnetization)
    3286            0 :          DO l = 0, MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
    3287            0 :             ll = 2*(2*l + 1)
    3288            0 :             ii = 0
    3289            0 :             DO i = 1, SIZE(ncalc, 2)
    3290            0 :                IF (ncalc(l, i) == 0) CYCLE
    3291            0 :                IF (ncalc(l, i) == ll) CYCLE
    3292            0 :                IF (ncalc(l, i) > dmag .AND. (ll - ncalc(l, i)) > dmag) THEN
    3293              :                   ii = i
    3294              :                   EXIT
    3295              :                END IF
    3296              :             END DO
    3297            0 :             IF (ii /= 0) THEN
    3298            0 :                edelta(l, ii, 1) = magnetization*0.5_dp
    3299            0 :                edelta(l, ii, 2) = -magnetization*0.5_dp
    3300            0 :                EXIT
    3301              :             END IF
    3302              :          END DO
    3303            0 :          IF (ii == 0) THEN
    3304              :             CALL cp_abort(__LOCATION__, &
    3305            0 :                           "Magnetization value cannot be imposed for this atom type")
    3306              :          END IF
    3307              :       END IF
    3308              : 
    3309        45297 :       IF (qs_kind%ghost .OR. qs_kind%floating) THEN
    3310          404 :          nelem = 0
    3311          404 :          ncore = 0
    3312          404 :          ncalc = 0
    3313          404 :          edelta = 0.0_dp
    3314              :       END IF
    3315              : 
    3316        45297 :    END SUBROUTINE init_atom_electronic_state
    3317              : 
    3318              : ! **************************************************************************************************
    3319              : !> \brief ...
    3320              : !> \param econf ...
    3321              : !> \param z ...
    3322              : !> \param ncalc ...
    3323              : !> \param ncore ...
    3324              : !> \param nelem ...
    3325              : ! **************************************************************************************************
    3326        24309 :    SUBROUTINE set_pseudo_state(econf, z, ncalc, ncore, nelem)
    3327              :       INTEGER, DIMENSION(:), POINTER                     :: econf
    3328              :       INTEGER, INTENT(IN)                                :: z
    3329              :       INTEGER, DIMENSION(0:lmat, 10), INTENT(OUT)        :: ncalc, ncore, nelem
    3330              : 
    3331              :       CHARACTER(LEN=default_string_length)               :: message
    3332              :       INTEGER                                            :: ii, iounit, l, ll, lmin, nc, nn
    3333              :       INTEGER, DIMENSION(0:lmat)                         :: econfx
    3334              :       TYPE(cp_logger_type), POINTER                      :: logger
    3335              : 
    3336        24309 :       NULLIFY (logger)
    3337        24309 :       logger => cp_get_default_logger()
    3338        24309 :       iounit = cp_logger_get_default_io_unit(logger)
    3339              : 
    3340        24309 :       econfx = 0
    3341        66584 :       econfx(0:SIZE(econf) - 1) = econf
    3342        66584 :       IF (SUM(econf) >= 0) THEN
    3343        66516 :          lmin = MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
    3344              :          ! number of core electrons
    3345        66516 :          nc = z - SUM(econf)
    3346              :          ! setup ncore
    3347        24275 :          ncore = 0
    3348         9062 :          SELECT CASE (nc)
    3349              :          CASE (0)
    3350              :          CASE (2)
    3351         9062 :             ncore(0, 1) = 2
    3352              :          CASE (10)
    3353         2232 :             ncore(0, 1) = 2
    3354         2232 :             ncore(0, 2) = 2
    3355         2232 :             ncore(1, 1) = 6
    3356              :          CASE (18)
    3357           58 :             ncore(0, 1) = 2
    3358           58 :             ncore(0, 2) = 2
    3359           58 :             ncore(0, 3) = 2
    3360           58 :             ncore(1, 1) = 6
    3361           58 :             ncore(1, 2) = 6
    3362              :          CASE (28)
    3363           16 :             ncore(0, 1) = 2
    3364           16 :             ncore(0, 2) = 2
    3365           16 :             ncore(0, 3) = 2
    3366           16 :             ncore(1, 1) = 6
    3367           16 :             ncore(1, 2) = 6
    3368           16 :             ncore(2, 1) = 10
    3369              :          CASE (36)
    3370            0 :             ncore(0, 1) = 2
    3371            0 :             ncore(0, 2) = 2
    3372            0 :             ncore(0, 3) = 2
    3373            0 :             ncore(0, 4) = 2
    3374            0 :             ncore(1, 1) = 6
    3375            0 :             ncore(1, 2) = 6
    3376            0 :             ncore(1, 3) = 6
    3377            0 :             ncore(2, 1) = 10
    3378              :          CASE (46)
    3379           36 :             ncore(0, 1) = 2
    3380           36 :             ncore(0, 2) = 2
    3381           36 :             ncore(0, 3) = 2
    3382           36 :             ncore(0, 4) = 2
    3383           36 :             ncore(1, 1) = 6
    3384           36 :             ncore(1, 2) = 6
    3385           36 :             ncore(1, 3) = 6
    3386           36 :             ncore(2, 1) = 10
    3387           36 :             ncore(2, 2) = 10
    3388              :          CASE (54)
    3389            4 :             ncore(0, 1) = 2
    3390            4 :             ncore(0, 2) = 2
    3391            4 :             ncore(0, 3) = 2
    3392            4 :             ncore(0, 4) = 2
    3393            4 :             ncore(0, 5) = 2
    3394            4 :             ncore(1, 1) = 6
    3395            4 :             ncore(1, 2) = 6
    3396            4 :             ncore(1, 3) = 6
    3397            4 :             ncore(1, 4) = 6
    3398            4 :             ncore(2, 1) = 10
    3399            4 :             ncore(2, 2) = 10
    3400              :          CASE (60)
    3401           30 :             ncore(0, 1) = 2
    3402           30 :             ncore(0, 2) = 2
    3403           30 :             ncore(0, 3) = 2
    3404           30 :             ncore(0, 4) = 2
    3405           30 :             ncore(1, 1) = 6
    3406           30 :             ncore(1, 2) = 6
    3407           30 :             ncore(1, 3) = 6
    3408           30 :             ncore(2, 1) = 10
    3409           30 :             ncore(2, 2) = 10
    3410           30 :             ncore(3, 1) = 14
    3411              :          CASE (68)
    3412          246 :             ncore(0, 1) = 2
    3413          246 :             ncore(0, 2) = 2
    3414          246 :             ncore(0, 3) = 2
    3415          246 :             ncore(0, 4) = 2
    3416          246 :             ncore(0, 5) = 2
    3417          246 :             ncore(1, 1) = 6
    3418          246 :             ncore(1, 2) = 6
    3419          246 :             ncore(1, 3) = 6
    3420          246 :             ncore(1, 4) = 6
    3421          246 :             ncore(2, 1) = 10
    3422          246 :             ncore(2, 2) = 10
    3423          246 :             ncore(3, 1) = 14
    3424              :          CASE (78)
    3425           12 :             ncore(0, 1) = 2
    3426           12 :             ncore(0, 2) = 2
    3427           12 :             ncore(0, 3) = 2
    3428           12 :             ncore(0, 4) = 2
    3429           12 :             ncore(0, 5) = 2
    3430           12 :             ncore(1, 1) = 6
    3431           12 :             ncore(1, 2) = 6
    3432           12 :             ncore(1, 3) = 6
    3433           12 :             ncore(1, 4) = 6
    3434           12 :             ncore(2, 1) = 10
    3435           12 :             ncore(2, 2) = 10
    3436           12 :             ncore(2, 3) = 10
    3437           12 :             ncore(3, 1) = 14
    3438              :             ! 79 - 92 5f incore PP
    3439              :          CASE (79)
    3440            0 :             ncore(0, 1) = 2
    3441            0 :             ncore(0, 2) = 2
    3442            0 :             ncore(0, 3) = 2
    3443            0 :             ncore(0, 4) = 2
    3444            0 :             ncore(0, 5) = 2
    3445            0 :             ncore(1, 1) = 6
    3446            0 :             ncore(1, 2) = 6
    3447            0 :             ncore(1, 3) = 6
    3448            0 :             ncore(1, 4) = 6
    3449            0 :             ncore(2, 1) = 10
    3450            0 :             ncore(2, 2) = 10
    3451            0 :             ncore(2, 3) = 10
    3452            0 :             ncore(3, 1) = 14
    3453            0 :             ncore(3, 2) = 1
    3454              :          CASE (80)
    3455            0 :             ncore(0, 1) = 2
    3456            0 :             ncore(0, 2) = 2
    3457            0 :             ncore(0, 3) = 2
    3458            0 :             ncore(0, 4) = 2
    3459            0 :             ncore(0, 5) = 2
    3460            0 :             ncore(1, 1) = 6
    3461            0 :             ncore(1, 2) = 6
    3462            0 :             ncore(1, 3) = 6
    3463            0 :             ncore(1, 4) = 6
    3464            0 :             ncore(2, 1) = 10
    3465            0 :             ncore(2, 2) = 10
    3466            0 :             ncore(2, 3) = 10
    3467            0 :             ncore(3, 1) = 14
    3468            0 :             ncore(3, 2) = 2
    3469              :          CASE (81)
    3470            0 :             ncore(0, 1) = 2
    3471            0 :             ncore(0, 2) = 2
    3472            0 :             ncore(0, 3) = 2
    3473            0 :             ncore(0, 4) = 2
    3474            0 :             ncore(0, 5) = 2
    3475            0 :             ncore(1, 1) = 6
    3476            0 :             ncore(1, 2) = 6
    3477            0 :             ncore(1, 3) = 6
    3478            0 :             ncore(1, 4) = 6
    3479            0 :             ncore(2, 1) = 10
    3480            0 :             ncore(2, 2) = 10
    3481            0 :             ncore(2, 3) = 10
    3482            0 :             ncore(3, 1) = 14
    3483            0 :             ncore(3, 2) = 3
    3484              :          CASE (82)
    3485            0 :             ncore(0, 1) = 2
    3486            0 :             ncore(0, 2) = 2
    3487            0 :             ncore(0, 3) = 2
    3488            0 :             ncore(0, 4) = 2
    3489            0 :             ncore(0, 5) = 2
    3490            0 :             ncore(1, 1) = 6
    3491            0 :             ncore(1, 2) = 6
    3492            0 :             ncore(1, 3) = 6
    3493            0 :             ncore(1, 4) = 6
    3494            0 :             ncore(2, 1) = 10
    3495            0 :             ncore(2, 2) = 10
    3496            0 :             ncore(2, 3) = 10
    3497            0 :             ncore(3, 1) = 14
    3498            0 :             ncore(3, 2) = 4
    3499              :          CASE (83)
    3500            0 :             ncore(0, 1) = 2
    3501            0 :             ncore(0, 2) = 2
    3502            0 :             ncore(0, 3) = 2
    3503            0 :             ncore(0, 4) = 2
    3504            0 :             ncore(0, 5) = 2
    3505            0 :             ncore(1, 1) = 6
    3506            0 :             ncore(1, 2) = 6
    3507            0 :             ncore(1, 3) = 6
    3508            0 :             ncore(1, 4) = 6
    3509            0 :             ncore(2, 1) = 10
    3510            0 :             ncore(2, 2) = 10
    3511            0 :             ncore(2, 3) = 10
    3512            0 :             ncore(3, 1) = 14
    3513            0 :             ncore(3, 2) = 5
    3514              :          CASE (84)
    3515            0 :             ncore(0, 1) = 2
    3516            0 :             ncore(0, 2) = 2
    3517            0 :             ncore(0, 3) = 2
    3518            0 :             ncore(0, 4) = 2
    3519            0 :             ncore(0, 5) = 2
    3520            0 :             ncore(1, 1) = 6
    3521            0 :             ncore(1, 2) = 6
    3522            0 :             ncore(1, 3) = 6
    3523            0 :             ncore(1, 4) = 6
    3524            0 :             ncore(2, 1) = 10
    3525            0 :             ncore(2, 2) = 10
    3526            0 :             ncore(2, 3) = 10
    3527            0 :             ncore(3, 1) = 14
    3528            0 :             ncore(3, 2) = 6
    3529              :          CASE (85)
    3530            0 :             ncore(0, 1) = 2
    3531            0 :             ncore(0, 2) = 2
    3532            0 :             ncore(0, 3) = 2
    3533            0 :             ncore(0, 4) = 2
    3534            0 :             ncore(0, 5) = 2
    3535            0 :             ncore(1, 1) = 6
    3536            0 :             ncore(1, 2) = 6
    3537            0 :             ncore(1, 3) = 6
    3538            0 :             ncore(1, 4) = 6
    3539            0 :             ncore(2, 1) = 10
    3540            0 :             ncore(2, 2) = 10
    3541            0 :             ncore(2, 3) = 10
    3542            0 :             ncore(3, 1) = 14
    3543            0 :             ncore(3, 2) = 7
    3544              :          CASE (86)
    3545              :             ! this is not Rn core, add double assignment below
    3546            0 :             ncore(0, 1) = 2
    3547            0 :             ncore(0, 2) = 2
    3548            0 :             ncore(0, 3) = 2
    3549            0 :             ncore(0, 4) = 2
    3550            0 :             ncore(0, 5) = 2
    3551            0 :             ncore(1, 1) = 6
    3552            0 :             ncore(1, 2) = 6
    3553            0 :             ncore(1, 3) = 6
    3554            0 :             ncore(1, 4) = 6
    3555            0 :             ncore(2, 1) = 10
    3556            0 :             ncore(2, 2) = 10
    3557            0 :             ncore(2, 3) = 10
    3558            0 :             ncore(3, 1) = 14
    3559            0 :             ncore(3, 2) = 8
    3560              :          CASE (87)
    3561            0 :             ncore(0, 1) = 2
    3562            0 :             ncore(0, 2) = 2
    3563            0 :             ncore(0, 3) = 2
    3564            0 :             ncore(0, 4) = 2
    3565            0 :             ncore(0, 5) = 2
    3566            0 :             ncore(1, 1) = 6
    3567            0 :             ncore(1, 2) = 6
    3568            0 :             ncore(1, 3) = 6
    3569            0 :             ncore(1, 4) = 6
    3570            0 :             ncore(2, 1) = 10
    3571            0 :             ncore(2, 2) = 10
    3572            0 :             ncore(2, 3) = 10
    3573            0 :             ncore(3, 1) = 14
    3574            0 :             ncore(3, 2) = 9
    3575              :          CASE (88)
    3576            0 :             ncore(0, 1) = 2
    3577            0 :             ncore(0, 2) = 2
    3578            0 :             ncore(0, 3) = 2
    3579            0 :             ncore(0, 4) = 2
    3580            0 :             ncore(0, 5) = 2
    3581            0 :             ncore(1, 1) = 6
    3582            0 :             ncore(1, 2) = 6
    3583            0 :             ncore(1, 3) = 6
    3584            0 :             ncore(1, 4) = 6
    3585            0 :             ncore(2, 1) = 10
    3586            0 :             ncore(2, 2) = 10
    3587            0 :             ncore(2, 3) = 10
    3588            0 :             ncore(3, 1) = 14
    3589            0 :             ncore(3, 2) = 10
    3590              :          CASE (89)
    3591            0 :             ncore(0, 1) = 2
    3592            0 :             ncore(0, 2) = 2
    3593            0 :             ncore(0, 3) = 2
    3594            0 :             ncore(0, 4) = 2
    3595            0 :             ncore(0, 5) = 2
    3596            0 :             ncore(1, 1) = 6
    3597            0 :             ncore(1, 2) = 6
    3598            0 :             ncore(1, 3) = 6
    3599            0 :             ncore(1, 4) = 6
    3600            0 :             ncore(2, 1) = 10
    3601            0 :             ncore(2, 2) = 10
    3602            0 :             ncore(2, 3) = 10
    3603            0 :             ncore(3, 1) = 14
    3604            0 :             ncore(3, 2) = 11
    3605              :          CASE (90)
    3606            0 :             ncore(0, 1) = 2
    3607            0 :             ncore(0, 2) = 2
    3608            0 :             ncore(0, 3) = 2
    3609            0 :             ncore(0, 4) = 2
    3610            0 :             ncore(0, 5) = 2
    3611            0 :             ncore(1, 1) = 6
    3612            0 :             ncore(1, 2) = 6
    3613            0 :             ncore(1, 3) = 6
    3614            0 :             ncore(1, 4) = 6
    3615            0 :             ncore(2, 1) = 10
    3616            0 :             ncore(2, 2) = 10
    3617            0 :             ncore(2, 3) = 10
    3618            0 :             ncore(3, 1) = 14
    3619            0 :             ncore(3, 2) = 12
    3620              :          CASE (91)
    3621            0 :             ncore(0, 1) = 2
    3622            0 :             ncore(0, 2) = 2
    3623            0 :             ncore(0, 3) = 2
    3624            0 :             ncore(0, 4) = 2
    3625            0 :             ncore(0, 5) = 2
    3626            0 :             ncore(1, 1) = 6
    3627            0 :             ncore(1, 2) = 6
    3628            0 :             ncore(1, 3) = 6
    3629            0 :             ncore(1, 4) = 6
    3630            0 :             ncore(2, 1) = 10
    3631            0 :             ncore(2, 2) = 10
    3632            0 :             ncore(2, 3) = 10
    3633            0 :             ncore(3, 1) = 14
    3634            0 :             ncore(3, 2) = 13
    3635              :          CASE (92)
    3636            0 :             ncore(0, 1) = 2
    3637            0 :             ncore(0, 2) = 2
    3638            0 :             ncore(0, 3) = 2
    3639            0 :             ncore(0, 4) = 2
    3640            0 :             ncore(0, 5) = 2
    3641            0 :             ncore(1, 1) = 6
    3642            0 :             ncore(1, 2) = 6
    3643            0 :             ncore(1, 3) = 6
    3644            0 :             ncore(1, 4) = 6
    3645            0 :             ncore(2, 1) = 10
    3646            0 :             ncore(2, 2) = 10
    3647            0 :             ncore(2, 3) = 10
    3648            0 :             ncore(3, 1) = 14
    3649            0 :             ncore(3, 2) = 14
    3650              :          CASE DEFAULT
    3651        24275 :             ncore(0, 1) = -1
    3652              :          END SELECT
    3653              :          ! special cases of double assignments
    3654        24275 :          IF (z == 65 .AND. econfx(3) == 0) THEN
    3655              :             ! 4f in core for Tb
    3656            4 :             ncore = 0
    3657            4 :             ncore(0, 1) = -1
    3658              :          END IF
    3659              :          ! if there is still no core, check for special cases
    3660        24275 :          IF (ncore(0, 1) <= 0) THEN
    3661        12583 :             IF (z >= 58 .AND. z <= 71) THEN
    3662              :                ! 4f-in-core PPs for lanthanides
    3663          280 :                nc = z - SUM(econf)
    3664              :                ! setup ncore
    3665           56 :                ncore = 0
    3666            0 :                SELECT CASE (nc)
    3667              :                CASE (29:42)
    3668            0 :                   ncore(0, 1) = 2
    3669            0 :                   ncore(0, 2) = 2
    3670            0 :                   ncore(0, 3) = 2
    3671            0 :                   ncore(1, 1) = 6
    3672            0 :                   ncore(1, 2) = 6
    3673            0 :                   ncore(2, 1) = 10
    3674            0 :                   ncore(3, 1) = nc - 28
    3675              :                   message = "A small-core pseudopotential with 4f-in-core is used for the lanthanide "// &
    3676            0 :                             TRIM(ptable(z)%symbol)
    3677            0 :                   CPHINT(TRIM(message))
    3678              :                CASE (47:60)
    3679           56 :                   ncore(0, 1) = 2
    3680           56 :                   ncore(0, 2) = 2
    3681           56 :                   ncore(0, 3) = 2
    3682           56 :                   ncore(0, 4) = 2
    3683           56 :                   ncore(1, 1) = 6
    3684           56 :                   ncore(1, 2) = 6
    3685           56 :                   ncore(1, 3) = 6
    3686           56 :                   ncore(2, 1) = 10
    3687           56 :                   ncore(2, 2) = 10
    3688           56 :                   ncore(3, 1) = nc - 46
    3689              :                   message = "A medium-core pseudopotential with 4f-in-core is used for the lanthanide "// &
    3690           56 :                             TRIM(ptable(z)%symbol)
    3691           56 :                   CPHINT(TRIM(message))
    3692              :                CASE DEFAULT
    3693           56 :                   ncore(0, 1) = -1
    3694              :                END SELECT
    3695              :             END IF
    3696              :          END IF
    3697              :          ! if the core is established, finish the setup
    3698        24275 :          IF (ncore(0, 1) >= 0) THEN
    3699       121375 :             DO l = 0, lmin
    3700        97100 :                ll = 2*(2*l + 1)
    3701      1068100 :                nn = SUM(ncore(l, :)) + econfx(l)
    3702        97100 :                ii = 0
    3703        24275 :                DO
    3704       117092 :                   ii = ii + 1
    3705       117092 :                   IF (nn <= ll) THEN
    3706        97100 :                      nelem(l, ii) = nn
    3707              :                      EXIT
    3708              :                   ELSE
    3709        19992 :                      nelem(l, ii) = ll
    3710        19992 :                      nn = nn - ll
    3711              :                   END IF
    3712              :                END DO
    3713              :             END DO
    3714      1723525 :             ncalc = nelem - ncore
    3715              :          ELSE
    3716              :             ! test for compatibility of valence occupation and full atomic occupation
    3717            0 :             IF (iounit > 0) THEN
    3718            0 :                WRITE (iounit, "(/,A,A2)") "WARNING: Core states irregular for atom type ", ptable(z)%symbol
    3719            0 :                WRITE (iounit, "(A,10I3)") "WARNING: Redefine ELEC_CONF in the KIND section"
    3720            0 :                CPABORT("Incompatible Atomic Occupations Detected")
    3721              :             END IF
    3722              :          END IF
    3723              :       ELSE
    3724           34 :          lmin = MIN(lmat, UBOUND(ptable(z)%e_conv, 1))
    3725           34 :          ncore = 0
    3726           34 :          ncalc = 0
    3727          170 :          DO l = 0, lmin
    3728          136 :             ll = 2*(2*l + 1)
    3729          136 :             nn = ABS(econfx(l))
    3730          136 :             ii = 0
    3731           34 :             DO
    3732          136 :                ii = ii + 1
    3733          136 :                IF (nn <= ll) THEN
    3734          136 :                   ncalc(l, ii) = -nn
    3735              :                   EXIT
    3736              :                ELSE
    3737            0 :                   ncalc(l, ii) = -ll
    3738            0 :                   nn = nn - ll
    3739              :                END IF
    3740              :             END DO
    3741              :          END DO
    3742           34 :          nelem = ncalc
    3743              :       END IF
    3744              : 
    3745        24309 :    END SUBROUTINE set_pseudo_state
    3746              : 
    3747              : ! **************************************************************************************************
    3748              : !> \brief finds if a given qs run needs to use nlcc
    3749              : !> \param qs_kind_set ...
    3750              : !> \return ...
    3751              : ! **************************************************************************************************
    3752        29334 :    FUNCTION has_nlcc(qs_kind_set) RESULT(nlcc)
    3753              : 
    3754              :       TYPE(qs_kind_type), DIMENSION(:)                   :: qs_kind_set
    3755              :       LOGICAL                                            :: nlcc
    3756              : 
    3757              :       INTEGER                                            :: ikind
    3758              :       LOGICAL                                            :: nlcc_present
    3759              :       TYPE(gth_potential_type), POINTER                  :: gth_potential
    3760              :       TYPE(sgp_potential_type), POINTER                  :: sgp_potential
    3761              : 
    3762        29334 :       nlcc = .FALSE.
    3763              : 
    3764        87739 :       DO ikind = 1, SIZE(qs_kind_set)
    3765        58405 :          CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential, sgp_potential=sgp_potential)
    3766        87739 :          IF (ASSOCIATED(gth_potential)) THEN
    3767        35479 :             CALL get_potential(potential=gth_potential, nlcc_present=nlcc_present)
    3768        35479 :             nlcc = nlcc .OR. nlcc_present
    3769        22926 :          ELSEIF (ASSOCIATED(sgp_potential)) THEN
    3770          336 :             CALL get_potential(potential=sgp_potential, has_nlcc=nlcc_present)
    3771          336 :             nlcc = nlcc .OR. nlcc_present
    3772              :          END IF
    3773              :       END DO
    3774              : 
    3775        29334 :    END FUNCTION has_nlcc
    3776              : 
    3777              : ! **************************************************************************************************
    3778              : 
    3779            0 : END MODULE qs_kind_types
        

Generated by: LCOV version 2.0-1