cp_control_utils.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief Utilities to set up the control types
! **************************************************************************************************
MODULE cp_control_utils
   USE bibliography,                    ONLY: &
        andreussi2012, dewar1977, dewar1985, elstner1998, fattebert2002, grimme2017, hu2007, &
        katbashev2025, krack2000, lippert1997, lippert1999, porezag1995, pracht2019, repasky2002, &
        rocha2006, schenter2008, seifert1996, souza2002, stengel2009, stewart1989, stewart2007, &
        thiel1992, umari2002, vanvoorhis2015, vandevondele2005a, vandevondele2005b, yin2017, &
        zhechkov2005, cite_reference
   USE cell_types,                      ONLY: cell_transform_input_cartesian,&
                                              cell_type
   USE cp_control_types,                ONLY: &
        admm_control_create, admm_control_type, ddapc_control_create, ddapc_restraint_type, &
        dft_control_create, dft_control_type, efield_type, expot_control_create, &
        maxwell_control_create, qs_control_type, rixs_control_type, tddfpt2_control_type, &
        xtb_control_type, xtb_reference_cli_type
   USE cp_files,                        ONLY: close_file,&
                                              open_file
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
                                              cp_print_key_unit_nr
   USE cp_parser_methods,               ONLY: parser_read_line
   USE cp_parser_types,                 ONLY: cp_parser_type,&
                                              parser_create,&
                                              parser_release,&
                                              parser_reset
   USE cp_units,                        ONLY: cp_unit_from_cp2k,&
                                              cp_unit_to_cp2k
   USE eeq_input,                       ONLY: read_eeq_param
   USE force_fields_input,              ONLY: read_gp_section
   USE input_constants,                 ONLY: &
        admm1_type, admm2_type, admmp_type, admmq_type, admms_type, constant_env, custom_env, &
        do_admm_aux_exch_func_bee, do_admm_aux_exch_func_bee_libxc, do_admm_aux_exch_func_default, &
        do_admm_aux_exch_func_default_libxc, do_admm_aux_exch_func_none, &
        do_admm_aux_exch_func_opt, do_admm_aux_exch_func_opt_libxc, do_admm_aux_exch_func_pbex, &
        do_admm_aux_exch_func_pbex_libxc, do_admm_aux_exch_func_sx_libxc, &
        do_admm_basis_projection, do_admm_blocked_projection, do_admm_blocking_purify_full, &
        do_admm_charge_constrained_projection, do_admm_exch_scaling_merlot, &
        do_admm_exch_scaling_none, do_admm_purify_cauchy, do_admm_purify_cauchy_subspace, &
        do_admm_purify_mcweeny, do_admm_purify_mo_diag, do_admm_purify_mo_no_diag, &
        do_admm_purify_none, do_admm_purify_none_dm, do_ddapc_constraint, do_ddapc_restraint, &
        do_method_am1, do_method_dftb, do_method_gapw, do_method_gapw_xc, do_method_gpw, &
        do_method_lrigpw, do_method_mndo, do_method_mndod, do_method_ofgpw, do_method_pdg, &
        do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rigpw, &
        do_method_rm1, do_method_xtb, do_pwgrid_ns_fullspace, do_pwgrid_ns_halfspace, &
        do_pwgrid_spherical, do_s2_constraint, do_s2_restraint, do_se_is_kdso, do_se_is_kdso_d, &
        do_se_is_slater, do_se_lr_ewald, do_se_lr_ewald_gks, do_se_lr_ewald_r3, do_se_lr_none, &
        gapw_1c_large, gapw_1c_medium, gapw_1c_orb, gapw_1c_small, gapw_1c_very_large, &
        gaussian_env, gfn1xtb, gfn_tblite, kg_tnadd_embed, kg_tnadd_embed_ri, no_admm_type, &
        numerical, ramp_env, real_time_propagation, rtp_method_bse, sccs_andreussi, &
        sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7, sccs_derivative_fft, &
        sccs_fattebert_gygi, sic_ad, sic_eo, sic_list_all, sic_list_unpaired, sic_mauri_spz, &
        sic_mauri_us, sic_none, slater, tblite_cli_born_kernel_auto, &
        tblite_cli_solution_state_gsolv, tblite_cli_solvation_alpb, tblite_cli_solvation_cpcm, &
        tblite_cli_solvation_gb, tblite_cli_solvation_gbe, tblite_cli_solvation_gbsa, &
        tblite_guess_ceh, tblite_mixer_memory_inherit, tblite_scc_mixer_auto, &
        tblite_scc_mixer_cp2k, tblite_scc_mixer_none, tblite_scc_mixer_tblite, tblite_solver_gvd, &
        tblite_solver_gvr, tddfpt_dipole_length, tddfpt_kernel_stda, use_mom_ref_user, &
        xtb_vdw_type_d3, xtb_vdw_type_d4, xtb_vdw_type_none
   USE input_cp2k_check,                ONLY: xc_functionals_expand
   USE input_cp2k_dft,                  ONLY: create_dft_section
   USE input_enumeration_types,         ONLY: enum_i2c,&
                                              enumeration_type
   USE input_keyword_types,             ONLY: keyword_get,&
                                              keyword_type
   USE input_section_types,             ONLY: &
        section_get_ival, section_get_keyword, section_release, section_type, section_vals_get, &
        section_vals_get_subs_vals, section_vals_type, section_vals_val_get, section_vals_val_set
   USE kinds,                           ONLY: default_path_length,&
                                              default_string_length,&
                                              dp
   USE mathconstants,                   ONLY: fourpi
   USE pair_potential_types,            ONLY: pair_potential_reallocate
   USE periodic_table,                  ONLY: get_ptable_info
   USE qs_cdft_utils,                   ONLY: read_cdft_control_section
   USE smeagol_control_types,           ONLY: read_smeagol_control
   USE string_utilities,                ONLY: uppercase
   USE util,                            ONLY: sort
   USE xas_tdp_types,                   ONLY: read_xas_tdp_control
   USE xc,                              ONLY: xc_uses_kinetic_energy_density,&
                                              xc_uses_norm_drho
   USE xc_input_constants,              ONLY: xc_deriv_collocate
   USE xc_write_output,                 ONLY: xc_write
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
 
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_control_utils'
 
   PUBLIC :: read_dft_control, &
             read_rixs_control, &
             read_mgrid_section, &
             read_qs_section, &
             read_tddfpt2_control, &
             write_dft_control, &
             write_qs_control, &
             write_admm_control, &
             read_ddapc_section
CONTAINS
 
! **************************************************************************************************
!> \brief ...
!> \param dft_control ...
!> \param dft_section ...
!> \param cell ...
! **************************************************************************************************
   SUBROUTINE read_dft_control(dft_control, dft_section, cell)
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(section_vals_type), POINTER                   :: dft_section
      TYPE(cell_type), OPTIONAL, POINTER                 :: cell
 
      CHARACTER(len=default_path_length)                 :: basis_set_file_name, &
                                                            intensities_file_name, &
                                                            potential_file_name
      CHARACTER(LEN=default_string_length), &
         DIMENSION(:), POINTER                           :: tmpstringlist
      INTEGER                                            :: admmtype, irep, isize, kg_tnadd_method, &
                                                            method_id, nrep, xc_deriv_method_id
      LOGICAL :: at_end, do_hfx, do_ot, do_rpa_admm, do_rtp, exopt1, exopt2, exopt3, explicit, &
         is_present, l_param, local_moment_possible, native_skala_grid, not_se, was_present
      REAL(kind=dp)                                      :: density_cut, gradient_cut, tau_cut
      REAL(kind=dp), DIMENSION(:), POINTER               :: pol
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(cp_parser_type)                               :: parser
      TYPE(section_vals_type), POINTER :: hairy_probes_section, hfx_section, kg_xc_fun_section, &
         kg_xc_section, maxwell_section, sccs_section, scf_section, tmp_section, xc_fun_section, &
         xc_gauxc_subsection, xc_section
 
      was_present = .false.
 
      logger => cp_get_default_logger()
 
      NULLIFY (kg_xc_fun_section, kg_xc_section, tmp_section, xc_fun_section, xc_section, xc_gauxc_subsection)
      ALLOCATE (dft_control)
      CALL dft_control_create(dft_control)
      ! determine wheather this is a semiempirical or DFTB run
      ! --> (no XC section needs to be provided)
      not_se = .true.
      CALL section_vals_val_get(dft_section, "QS%METHOD", i_val=method_id)
      SELECT CASE (method_id)
      CASE (do_method_dftb, do_method_xtb, do_method_mndo, do_method_am1, do_method_pm3, do_method_pnnl, &
            do_method_pm6, do_method_pm6fm, do_method_pdg, do_method_rm1, do_method_mndod)
         not_se = .false.
      END SELECT
      ! Check for XC section and XC_FUNCTIONAL section
      xc_section => section_vals_get_subs_vals(dft_section, "XC")
      CALL section_vals_get(xc_section, explicit=is_present)
      IF (.NOT. is_present .AND. not_se) THEN
         cpabort("XC section missing.")
      END IF
      IF (is_present) THEN
         CALL section_vals_val_get(xc_section, "density_cutoff", r_val=density_cut)
         CALL section_vals_val_get(xc_section, "gradient_cutoff", r_val=gradient_cut)
         CALL section_vals_val_get(xc_section, "tau_cutoff", r_val=tau_cut)
         ! Perform numerical stability checks and possibly correct the issues
         IF (density_cut <= epsilon(0.0_dp)*100.0_dp) &
            CALL cp_warn(__location__, &
                         "DENSITY_CUTOFF lower than 100*EPSILON, where EPSILON is the machine precision. "// &
                         "This may lead to numerical problems. Setting up shake_tol to 100*EPSILON! ")
         density_cut = max(epsilon(0.0_dp)*100.0_dp, density_cut)
         IF (gradient_cut <= epsilon(0.0_dp)*100.0_dp) &
            CALL cp_warn(__location__, &
                         "GRADIENT_CUTOFF lower than 100*EPSILON, where EPSILON is the machine precision. "// &
                         "This may lead to numerical problems. Setting up shake_tol to 100*EPSILON! ")
         gradient_cut = max(epsilon(0.0_dp)*100.0_dp, gradient_cut)
         IF (tau_cut <= epsilon(0.0_dp)*100.0_dp) &
            CALL cp_warn(__location__, &
                         "TAU_CUTOFF lower than 100*EPSILON, where EPSILON is the machine precision. "// &
                         "This may lead to numerical problems. Setting up shake_tol to 100*EPSILON! ")
         tau_cut = max(epsilon(0.0_dp)*100.0_dp, tau_cut)
         CALL section_vals_val_set(xc_section, "density_cutoff", r_val=density_cut)
         CALL section_vals_val_set(xc_section, "gradient_cutoff", r_val=gradient_cut)
         CALL section_vals_val_set(xc_section, "tau_cutoff", r_val=tau_cut)
      END IF
      xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
      CALL section_vals_get(xc_fun_section, explicit=is_present)
      IF (.NOT. is_present .AND. not_se) THEN
         cpabort("XC_FUNCTIONAL section missing.")
      END IF
 
      dft_control%use_gauxc = .false.
      IF (is_present) THEN
         xc_gauxc_subsection => section_vals_get_subs_vals(xc_fun_section, "GAUXC")
         CALL section_vals_get(xc_gauxc_subsection, explicit=dft_control%use_gauxc)
      END IF
 
      scf_section => section_vals_get_subs_vals(dft_section, "SCF")
      CALL section_vals_val_get(dft_section, "UKS", l_val=dft_control%uks)
      CALL section_vals_val_get(dft_section, "ROKS", l_val=dft_control%roks)
      IF (dft_control%uks .OR. dft_control%roks) THEN
         dft_control%nspins = 2
      ELSE
         dft_control%nspins = 1
      END IF
 
      dft_control%lsd = (dft_control%nspins > 1)
      dft_control%use_kinetic_energy_density = xc_uses_kinetic_energy_density(xc_fun_section, dft_control%lsd)
      IF (dft_control%use_gauxc) THEN
         native_skala_grid = .false.
         CALL section_vals_val_get(xc_gauxc_subsection, "NATIVE_GRID", l_val=native_skala_grid)
         IF (native_skala_grid) dft_control%use_kinetic_energy_density = .true.
      END IF
      tmp_section => section_vals_get_subs_vals(dft_section, "KG_METHOD")
      CALL section_vals_get(tmp_section, explicit=explicit)
      IF (explicit) THEN
         CALL section_vals_val_get(tmp_section, "TNADD_METHOD", i_val=kg_tnadd_method)
         IF (kg_tnadd_method == kg_tnadd_embed .OR. kg_tnadd_method == kg_tnadd_embed_ri) THEN
            kg_xc_section => section_vals_get_subs_vals(tmp_section, "XC")
            kg_xc_fun_section => section_vals_get_subs_vals(kg_xc_section, "XC_FUNCTIONAL")
            CALL section_vals_get(kg_xc_fun_section, explicit=is_present)
            IF (is_present) THEN
               dft_control%use_kinetic_energy_density = dft_control%use_kinetic_energy_density .OR. &
                                                        xc_uses_kinetic_energy_density(kg_xc_fun_section, &
                                                                                       dft_control%lsd)
            END IF
         END IF
      END IF
 
      xc_deriv_method_id = section_get_ival(xc_section, "XC_GRID%XC_DERIV")
      dft_control%drho_by_collocation = (xc_uses_norm_drho(xc_fun_section, dft_control%lsd) &
                                         .AND. (xc_deriv_method_id == xc_deriv_collocate))
      IF (dft_control%drho_by_collocation) THEN
         cpabort("derivatives by collocation not implemented")
      END IF
 
      ! Automatic auxiliary basis set generation
      CALL section_vals_val_get(dft_section, "AUTO_BASIS", n_rep_val=nrep)
      DO irep = 1, nrep
         CALL section_vals_val_get(dft_section, "AUTO_BASIS", i_rep_val=irep, c_vals=tmpstringlist)
         IF (SIZE(tmpstringlist) == 2) THEN
            CALL uppercase(tmpstringlist(2))
            SELECT CASE (tmpstringlist(2))
            CASE ("X")
               isize = -1
            CASE ("SMALL")
               isize = 0
            CASE ("MEDIUM")
               isize = 1
            CASE ("LARGE")
               isize = 2
            CASE ("HUGE")
               isize = 3
            CASE DEFAULT
               cpwarn("Unknown basis size in AUTO_BASIS keyword:"//trim(tmpstringlist(1)))
            END SELECT
            !
            SELECT CASE (tmpstringlist(1))
            CASE ("X")
            CASE ("RI_AUX")
               dft_control%auto_basis_ri_aux = isize
            CASE ("AUX_FIT")
               dft_control%auto_basis_aux_fit = isize
            CASE ("LRI_AUX")
               dft_control%auto_basis_lri_aux = isize
            CASE ("P_LRI_AUX")
               dft_control%auto_basis_p_lri_aux = isize
            CASE ("RI_HXC")
               dft_control%auto_basis_ri_hxc = isize
            CASE ("RI_XAS")
               dft_control%auto_basis_ri_xas = isize
            CASE ("RI_HFX")
               dft_control%auto_basis_ri_hfx = isize
            CASE DEFAULT
               cpwarn("Unknown basis type in AUTO_BASIS keyword:"//trim(tmpstringlist(1)))
            END SELECT
         ELSE
            CALL cp_abort(__location__, &
                          "AUTO_BASIS keyword in &DFT section has a wrong number of arguments.")
         END IF
      END DO
 
      !! check if we do wavefunction fitting
      tmp_section => section_vals_get_subs_vals(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD")
      CALL section_vals_get(tmp_section, explicit=is_present)
      !
      hfx_section => section_vals_get_subs_vals(xc_section, "HF")
      CALL section_vals_get(hfx_section, explicit=do_hfx)
      CALL section_vals_val_get(xc_section, "WF_CORRELATION%RI_RPA%ADMM", l_val=do_rpa_admm)
      is_present = is_present .AND. (do_hfx .OR. do_rpa_admm)
      !
      dft_control%do_admm = is_present
      dft_control%do_admm_mo = .false.
      dft_control%do_admm_dm = .false.
      IF (is_present) THEN
         do_ot = .false.
         CALL section_vals_val_get(scf_section, "OT%_SECTION_PARAMETERS_", l_val=do_ot)
         CALL admm_control_create(dft_control%admm_control)
 
         CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%ADMM_TYPE", i_val=admmtype)
         CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%ADMM_PURIFICATION_METHOD", explicit=exopt1)
         CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%METHOD", explicit=exopt2)
         CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%EXCH_SCALING_MODEL", explicit=exopt3)
         dft_control%admm_control%admm_type = admmtype
         SELECT CASE (admmtype)
         CASE (no_admm_type)
            CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%ADMM_PURIFICATION_METHOD", i_val=method_id)
            dft_control%admm_control%purification_method = method_id
            CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%METHOD", i_val=method_id)
            dft_control%admm_control%method = method_id
            CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%EXCH_SCALING_MODEL", i_val=method_id)
            dft_control%admm_control%scaling_model = method_id
         CASE (admm1_type)
            ! METHOD BASIS_PROJECTION
            ! ADMM_PURIFICATION_METHOD choose
            ! EXCH_SCALING_MODEL NONE
            CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%ADMM_PURIFICATION_METHOD", i_val=method_id)
            dft_control%admm_control%purification_method = method_id
            dft_control%admm_control%method = do_admm_basis_projection
            dft_control%admm_control%scaling_model = do_admm_exch_scaling_none
         CASE (admm2_type)
            ! METHOD BASIS_PROJECTION
            ! ADMM_PURIFICATION_METHOD NONE
            ! EXCH_SCALING_MODEL NONE
            dft_control%admm_control%purification_method = do_admm_purify_none
            dft_control%admm_control%method = do_admm_basis_projection
            dft_control%admm_control%scaling_model = do_admm_exch_scaling_none
         CASE (admms_type)
            ! ADMM_PURIFICATION_METHOD NONE
            ! METHOD CHARGE_CONSTRAINED_PROJECTION
            ! EXCH_SCALING_MODEL MERLOT
            dft_control%admm_control%purification_method = do_admm_purify_none
            dft_control%admm_control%method = do_admm_charge_constrained_projection
            dft_control%admm_control%scaling_model = do_admm_exch_scaling_merlot
         CASE (admmp_type)
            ! ADMM_PURIFICATION_METHOD NONE
            ! METHOD BASIS_PROJECTION
            ! EXCH_SCALING_MODEL MERLOT
            dft_control%admm_control%purification_method = do_admm_purify_none
            dft_control%admm_control%method = do_admm_basis_projection
            dft_control%admm_control%scaling_model = do_admm_exch_scaling_merlot
         CASE (admmq_type)
            ! ADMM_PURIFICATION_METHOD NONE
            ! METHOD CHARGE_CONSTRAINED_PROJECTION
            ! EXCH_SCALING_MODEL NONE
            dft_control%admm_control%purification_method = do_admm_purify_none
            dft_control%admm_control%method = do_admm_charge_constrained_projection
            dft_control%admm_control%scaling_model = do_admm_exch_scaling_none
         CASE DEFAULT
            CALL cp_abort(__location__, &
                          "ADMM_TYPE keyword in &AUXILIARY_DENSITY_MATRIX_METHOD section has a wrong value.")
         END SELECT
 
         CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%EPS_FILTER", &
                                   r_val=dft_control%admm_control%eps_filter)
 
         CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%EXCH_CORRECTION_FUNC", i_val=method_id)
         dft_control%admm_control%aux_exch_func = method_id
 
         ! parameters for X functional
         dft_control%admm_control%aux_exch_func_param = .false.
         CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%OPTX_A1", explicit=explicit, &
                                   r_val=dft_control%admm_control%aux_x_param(1))
         IF (explicit) dft_control%admm_control%aux_exch_func_param = .true.
         CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%OPTX_A2", explicit=explicit, &
                                   r_val=dft_control%admm_control%aux_x_param(2))
         IF (explicit) dft_control%admm_control%aux_exch_func_param = .true.
         CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%OPTX_GAMMA", explicit=explicit, &
                                   r_val=dft_control%admm_control%aux_x_param(3))
         IF (explicit) dft_control%admm_control%aux_exch_func_param = .true.
 
         CALL read_admm_block_list(dft_control%admm_control, dft_section)
 
         ! check for double assignments
         SELECT CASE (admmtype)
         CASE (admm2_type)
            IF (exopt2) CALL cp_warn(__location__, &
                                     "Value of ADMM_PURIFICATION_METHOD keyword will be overwritten with ADMM_TYPE selections.")
            IF (exopt3) CALL cp_warn(__location__, &
                                     "Value of EXCH_SCALING_MODEL keyword will be overwritten with ADMM_TYPE selections.")
         CASE (admm1_type, admms_type, admmp_type, admmq_type)
            IF (exopt1) CALL cp_warn(__location__, &
                                     "Value of METHOD keyword will be overwritten with ADMM_TYPE selections.")
            IF (exopt2) CALL cp_warn(__location__, &
                                     "Value of METHOD keyword will be overwritten with ADMM_TYPE selections.")
            IF (exopt3) CALL cp_warn(__location__, &
                                     "Value of EXCH_SCALING_MODEL keyword will be overwritten with ADMM_TYPE selections.")
         END SELECT
 
         !    In the case of charge-constrained projection (e.g. according to Merlot),
         !    there is no purification needed and hence, do_admm_purify_none has to be set.
 
         IF ((dft_control%admm_control%method == do_admm_blocking_purify_full .OR. &
              dft_control%admm_control%method == do_admm_blocked_projection) &
             .AND. dft_control%admm_control%scaling_model == do_admm_exch_scaling_merlot) THEN
            cpabort("ADMM: Blocking and Merlot scaling are mutually exclusive.")
         END IF
 
         IF (dft_control%admm_control%method == do_admm_charge_constrained_projection .AND. &
             dft_control%admm_control%purification_method /= do_admm_purify_none) THEN
            CALL cp_abort(__location__, &
                          "ADMM: In the case of METHOD=CHARGE_CONSTRAINED_PROJECTION, "// &
                          "ADMM_PURIFICATION_METHOD=NONE has to be set.")
         END IF
 
         IF (dft_control%admm_control%purification_method == do_admm_purify_mo_diag .OR. &
             dft_control%admm_control%purification_method == do_admm_purify_mo_no_diag) THEN
            IF (dft_control%admm_control%method /= do_admm_basis_projection) &
               cpabort("ADMM: Chosen purification requires BASIS_PROJECTION")
 
            IF (.NOT. do_ot) cpabort("ADMM: MO-based purification requires OT.")
         END IF
 
         IF (dft_control%admm_control%purification_method == do_admm_purify_none_dm .OR. &
             dft_control%admm_control%purification_method == do_admm_purify_mcweeny) THEN
            dft_control%do_admm_dm = .true.
         ELSE
            dft_control%do_admm_mo = .true.
         END IF
      END IF
 
      ! Set restricted to true, if both OT and ROKS are requested
      !MK in principle dft_control%restricted could be dropped completely like the
      !MK input key by using only dft_control%roks now
      CALL section_vals_val_get(scf_section, "OT%_SECTION_PARAMETERS_", l_val=l_param)
      dft_control%restricted = (dft_control%roks .AND. l_param)
 
      CALL section_vals_val_get(dft_section, "CHARGE", i_val=dft_control%charge)
      CALL section_vals_val_get(dft_section, "MULTIPLICITY", i_val=dft_control%multiplicity)
      CALL section_vals_val_get(dft_section, "RELAX_MULTIPLICITY", r_val=dft_control%relax_multiplicity)
      IF (dft_control%relax_multiplicity > 0.0_dp) THEN
         IF (.NOT. dft_control%uks) &
            CALL cp_abort(__location__, "The option RELAX_MULTIPLICITY is only valid for "// &
                          "unrestricted Kohn-Sham (UKS) calculations")
      END IF
 
      !Read the HAIR PROBES input section if present
      hairy_probes_section => section_vals_get_subs_vals(dft_section, "HAIRY_PROBES")
      CALL section_vals_get(hairy_probes_section, n_repetition=nrep, explicit=is_present)
 
      IF (is_present) THEN
         dft_control%hairy_probes = .true.
         ALLOCATE (dft_control%probe(nrep))
         CALL read_hairy_probes_sections(dft_control, hairy_probes_section)
      END IF
 
      ! check for the presence of the low spin roks section
      tmp_section => section_vals_get_subs_vals(dft_section, "LOW_SPIN_ROKS")
      CALL section_vals_get(tmp_section, explicit=dft_control%low_spin_roks)
 
      dft_control%sic_method_id = sic_none
      dft_control%sic_scaling_a = 1.0_dp
      dft_control%sic_scaling_b = 1.0_dp
 
      ! DFT+U
      dft_control%dft_plus_u = .false.
      CALL section_vals_val_get(dft_section, "PLUS_U_METHOD", i_val=method_id)
      dft_control%plus_u_method_id = method_id
 
      ! Smearing in use
      dft_control%smear = .false.
 
      ! Surface dipole correction
      dft_control%correct_surf_dip = .false.
      CALL section_vals_val_get(dft_section, "SURFACE_DIPOLE_CORRECTION", l_val=dft_control%correct_surf_dip)
      CALL section_vals_val_get(dft_section, "SURF_DIP_DIR", i_val=dft_control%dir_surf_dip)
      dft_control%pos_dir_surf_dip = -1.0_dp
      CALL section_vals_val_get(dft_section, "SURF_DIP_POS", r_val=dft_control%pos_dir_surf_dip)
      ! another logical variable, surf_dip_correct_switch, is introduced for
      ! implementation of "SURF_DIP_SWITCH" [SGh]
      dft_control%switch_surf_dip = .false.
      dft_control%surf_dip_correct_switch = dft_control%correct_surf_dip
      CALL section_vals_val_get(dft_section, "SURF_DIP_SWITCH", l_val=dft_control%switch_surf_dip)
      dft_control%correct_el_density_dip = .false.
      CALL section_vals_val_get(dft_section, "CORE_CORR_DIP", l_val=dft_control%correct_el_density_dip)
      IF (dft_control%correct_el_density_dip) THEN
         IF (dft_control%correct_surf_dip) THEN
            ! Do nothing, move on
         ELSE
            dft_control%correct_el_density_dip = .false.
            cpwarn("CORE_CORR_DIP keyword is activated only if SURFACE_DIPOLE_CORRECTION is TRUE")
         END IF
      END IF
 
      CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", &
                                c_val=basis_set_file_name)
      CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", &
                                c_val=potential_file_name)
 
      ! Read the input section
      tmp_section => section_vals_get_subs_vals(dft_section, "sic")
      CALL section_vals_val_get(tmp_section, "SIC_METHOD", &
                                i_val=dft_control%sic_method_id)
      CALL section_vals_val_get(tmp_section, "ORBITAL_SET", &
                                i_val=dft_control%sic_list_id)
      CALL section_vals_val_get(tmp_section, "SIC_SCALING_A", &
                                r_val=dft_control%sic_scaling_a)
      CALL section_vals_val_get(tmp_section, "SIC_SCALING_B", &
                                r_val=dft_control%sic_scaling_b)
 
      do_rtp = .false.
      tmp_section => section_vals_get_subs_vals(dft_section, "REAL_TIME_PROPAGATION")
      CALL section_vals_get(tmp_section, explicit=is_present)
      IF (is_present) THEN
         CALL read_rtp_section(dft_control, tmp_section)
         do_rtp = .true.
      END IF
 
      ! Read the input section
      tmp_section => section_vals_get_subs_vals(dft_section, "XAS")
      CALL section_vals_get(tmp_section, explicit=dft_control%do_xas_calculation)
      IF (dft_control%do_xas_calculation) THEN
         ! Override with section parameter
         CALL section_vals_val_get(tmp_section, "_SECTION_PARAMETERS_", &
                                   l_val=dft_control%do_xas_calculation)
      END IF
 
      tmp_section => section_vals_get_subs_vals(dft_section, "XAS_TDP")
      CALL section_vals_get(tmp_section, explicit=dft_control%do_xas_tdp_calculation)
      IF (dft_control%do_xas_tdp_calculation) THEN
         ! Override with section parameter
         CALL section_vals_val_get(tmp_section, "_SECTION_PARAMETERS_", &
                                   l_val=dft_control%do_xas_tdp_calculation)
      END IF
 
      ! Read the finite field input section
      dft_control%apply_efield = .false.
      dft_control%apply_efield_field = .false. !this is for RTP
      dft_control%apply_vector_potential = .false. !this is for RTP
      tmp_section => section_vals_get_subs_vals(dft_section, "EFIELD")
      CALL section_vals_get(tmp_section, n_repetition=nrep, explicit=is_present)
      IF (is_present) THEN
         ALLOCATE (dft_control%efield_fields(nrep))
         CALL read_efield_sections(dft_control, tmp_section, cell)
         IF (do_rtp) THEN
            IF (.NOT. dft_control%rtp_control%velocity_gauge) THEN
               dft_control%apply_efield_field = .true.
            ELSE
               dft_control%apply_vector_potential = .true.
               ! Use this input value of vector potential to (re)start RTP
               dft_control%rtp_control%vec_pot = dft_control%efield_fields(1)%efield%vec_pot_initial
            END IF
         ELSE
            dft_control%apply_efield = .true.
            cpassert(nrep == 1)
         END IF
      END IF
 
      ! Now, can try to guess polarisation in rtp
      IF (do_rtp) THEN
         ! tmp_section => section_vals_get_subs_vals(dft_section, "REAL_TIME_PROPAGATION%PRINT%POLARIZABILITY")
         ! CALL section_vals_get(tmp_section, explicit=is_present)
         local_moment_possible = (dft_control%rtp_control%rtp_method == rtp_method_bse) .OR. &
                                 ((.NOT. dft_control%rtp_control%periodic) .AND. dft_control%rtp_control%linear_scaling)
         IF (local_moment_possible .AND. (.NOT. ASSOCIATED(dft_control%rtp_control%print_pol_elements))) THEN
            tmp_section => section_vals_get_subs_vals(dft_section, "REAL_TIME_PROPAGATION")
            CALL guess_pol_elements(dft_control, &
                                    dft_control%rtp_control%print_pol_elements)
         END IF
      END IF
 
      ! Read the finite field input section for periodic fields
      tmp_section => section_vals_get_subs_vals(dft_section, "PERIODIC_EFIELD")
      CALL section_vals_get(tmp_section, explicit=dft_control%apply_period_efield)
      IF (dft_control%apply_period_efield) THEN
         ALLOCATE (dft_control%period_efield)
         CALL section_vals_val_get(tmp_section, "POLARISATION", r_vals=pol)
         dft_control%period_efield%polarisation(1:3) = pol(1:3)
         IF (PRESENT(cell)) THEN
            IF (ASSOCIATED(cell)) &
               CALL cell_transform_input_cartesian(cell, dft_control%period_efield%polarisation(1:3))
         END IF
         CALL section_vals_val_get(tmp_section, "D_FILTER", r_vals=pol)
         dft_control%period_efield%d_filter(1:3) = pol(1:3)
         IF (PRESENT(cell)) THEN
            IF (ASSOCIATED(cell)) CALL cell_transform_input_cartesian(cell, dft_control%period_efield%d_filter(1:3))
         END IF
         CALL section_vals_val_get(tmp_section, "INTENSITY", &
                                   r_val=dft_control%period_efield%strength)
         dft_control%period_efield%displacement_field = .false.
         CALL section_vals_val_get(tmp_section, "DISPLACEMENT_FIELD", &
                                   l_val=dft_control%period_efield%displacement_field)
 
         CALL section_vals_val_get(tmp_section, "INTENSITY_LIST", r_vals=pol)
 
         CALL section_vals_val_get(tmp_section, "INTENSITIES_FILE_NAME", c_val=intensities_file_name)
 
         IF (SIZE(pol) > 1 .OR. pol(1) /= 0.0_dp) THEN
            ! if INTENSITY_LIST is present, INTENSITY and INTENSITIES_FILE_NAME must not be present
            IF (dft_control%period_efield%strength /= 0.0_dp .OR. intensities_file_name /= "") THEN
               CALL cp_abort(__location__, "[PERIODIC FIELD] Only one of INTENSITY, INTENSITY_LIST "// &
                             "or INTENSITIES_FILE_NAME can be specified.")
            END IF
 
            ALLOCATE (dft_control%period_efield%strength_list(SIZE(pol)))
            dft_control%period_efield%strength_list(1:SIZE(pol)) = pol(1:SIZE(pol))
         END IF
 
         IF (intensities_file_name /= "") THEN
            ! if INTENSITIES_FILE_NAME is present, INTENSITY must not be present
            IF (dft_control%period_efield%strength /= 0.0_dp) THEN
               CALL cp_abort(__location__, "[PERIODIC FIELD] Only one of INTENSITY, INTENSITY_LIST "// &
                             "or INTENSITIES_FILE_NAME can be specified.")
            END IF
 
            CALL parser_create(parser, intensities_file_name)
 
            nrep = 0
            DO WHILE (.true.)
               CALL parser_read_line(parser, 1, at_end)
               IF (at_end) EXIT
               nrep = nrep + 1
            END DO
 
            IF (nrep == 0) THEN
               cpabort("[PERIODIC FIELD] No intensities found in INTENSITIES_FILE_NAME")
            END IF
 
            ALLOCATE (dft_control%period_efield%strength_list(nrep))
 
            CALL parser_reset(parser)
            DO irep = 1, nrep
               CALL parser_read_line(parser, 1)
               READ (parser%input_line, *) dft_control%period_efield%strength_list(irep)
            END DO
 
            CALL parser_release(parser)
         END IF
 
         CALL section_vals_val_get(tmp_section, "START_FRAME", &
                                   i_val=dft_control%period_efield%start_frame)
         CALL section_vals_val_get(tmp_section, "END_FRAME", &
                                   i_val=dft_control%period_efield%end_frame)
 
         IF (dft_control%period_efield%end_frame /= -1) THEN
            ! check if valid bounds are given
            ! if an end frame is given, the number of active frames must be a
            ! multiple of the number of intensities
            IF (dft_control%period_efield%start_frame > dft_control%period_efield%end_frame) THEN
               cpabort("[PERIODIC FIELD] START_FRAME > END_FRAME")
            ELSE IF (dft_control%period_efield%start_frame < 1) THEN
               cpabort("[PERIODIC FIELD] START_FRAME < 1")
            ELSE IF (mod(dft_control%period_efield%end_frame - &
                         dft_control%period_efield%start_frame + 1, SIZE(pol)) /= 0) THEN
               CALL cp_abort(__location__, &
                             "[PERIODIC FIELD] Number of active frames must be a multiple of the number of intensities")
            END IF
         END IF
 
         ! periodic fields don't work with RTP
         cpassert(.NOT. do_rtp)
         IF (dft_control%period_efield%displacement_field) THEN
            CALL cite_reference(stengel2009)
         ELSE
            CALL cite_reference(souza2002)
            CALL cite_reference(umari2002)
         END IF
      END IF
 
      ! Read the external potential input section
      tmp_section => section_vals_get_subs_vals(dft_section, "EXTERNAL_POTENTIAL")
      CALL section_vals_get(tmp_section, explicit=dft_control%apply_external_potential)
      IF (dft_control%apply_external_potential) THEN
         CALL expot_control_create(dft_control%expot_control)
         CALL section_vals_val_get(tmp_section, "READ_FROM_CUBE", &
                                   l_val=dft_control%expot_control%read_from_cube)
         CALL section_vals_val_get(tmp_section, "STATIC", &
                                   l_val=dft_control%expot_control%static)
         CALL section_vals_val_get(tmp_section, "SCALING_FACTOR", &
                                   r_val=dft_control%expot_control%scaling_factor)
         ! External potential using Maxwell equation
         maxwell_section => section_vals_get_subs_vals(tmp_section, "MAXWELL")
         CALL section_vals_get(maxwell_section, explicit=is_present)
         IF (is_present) THEN
            dft_control%expot_control%maxwell_solver = .true.
            CALL maxwell_control_create(dft_control%maxwell_control)
            ! read the input values from Maxwell section
            CALL section_vals_val_get(maxwell_section, "TEST_REAL", &
                                      r_val=dft_control%maxwell_control%real_test)
            CALL section_vals_val_get(maxwell_section, "TEST_INTEGER", &
                                      i_val=dft_control%maxwell_control%int_test)
            CALL section_vals_val_get(maxwell_section, "TEST_LOGICAL", &
                                      l_val=dft_control%maxwell_control%log_test)
         ELSE
            dft_control%expot_control%maxwell_solver = .false.
         END IF
      END IF
 
      ! Read the SCCS input section if present
      sccs_section => section_vals_get_subs_vals(dft_section, "SCCS")
      CALL section_vals_get(sccs_section, explicit=is_present)
      IF (is_present) THEN
         ! Check section parameter if SCCS is activated
         CALL section_vals_val_get(sccs_section, "_SECTION_PARAMETERS_", &
                                   l_val=dft_control%do_sccs)
         IF (dft_control%do_sccs) THEN
            ALLOCATE (dft_control%sccs_control)
            CALL section_vals_val_get(sccs_section, "RELATIVE_PERMITTIVITY", &
                                      r_val=dft_control%sccs_control%epsilon_solvent)
            CALL section_vals_val_get(sccs_section, "ALPHA", &
                                      r_val=dft_control%sccs_control%alpha_solvent)
            CALL section_vals_val_get(sccs_section, "BETA", &
                                      r_val=dft_control%sccs_control%beta_solvent)
            CALL section_vals_val_get(sccs_section, "DELTA_RHO", &
                                      r_val=dft_control%sccs_control%delta_rho)
            CALL section_vals_val_get(sccs_section, "DERIVATIVE_METHOD", &
                                      i_val=dft_control%sccs_control%derivative_method)
            CALL section_vals_val_get(sccs_section, "METHOD", &
                                      i_val=dft_control%sccs_control%method_id)
            CALL section_vals_val_get(sccs_section, "EPS_SCCS", &
                                      r_val=dft_control%sccs_control%eps_sccs)
            CALL section_vals_val_get(sccs_section, "EPS_SCF", &
                                      r_val=dft_control%sccs_control%eps_scf)
            CALL section_vals_val_get(sccs_section, "GAMMA", &
                                      r_val=dft_control%sccs_control%gamma_solvent)
            CALL section_vals_val_get(sccs_section, "MAX_ITER", &
                                      i_val=dft_control%sccs_control%max_iter)
            CALL section_vals_val_get(sccs_section, "MIXING", &
                                      r_val=dft_control%sccs_control%mixing)
            SELECT CASE (dft_control%sccs_control%method_id)
            CASE (sccs_andreussi)
               tmp_section => section_vals_get_subs_vals(sccs_section, "ANDREUSSI")
               CALL section_vals_val_get(tmp_section, "RHO_MAX", &
                                         r_val=dft_control%sccs_control%rho_max)
               CALL section_vals_val_get(tmp_section, "RHO_MIN", &
                                         r_val=dft_control%sccs_control%rho_min)
               IF (dft_control%sccs_control%rho_max < dft_control%sccs_control%rho_min) THEN
                  CALL cp_abort(__location__, &
                                "The SCCS parameter RHO_MAX is smaller than RHO_MIN. "// &
                                "Please, check your input!")
               END IF
               CALL cite_reference(andreussi2012)
            CASE (sccs_fattebert_gygi)
               tmp_section => section_vals_get_subs_vals(sccs_section, "FATTEBERT-GYGI")
               CALL section_vals_val_get(tmp_section, "BETA", &
                                         r_val=dft_control%sccs_control%beta)
               IF (dft_control%sccs_control%beta < 0.5_dp) THEN
                  CALL cp_abort(__location__, &
                                "A value smaller than 0.5 for the SCCS parameter beta "// &
                                "causes numerical problems. Please, check your input!")
               END IF
               CALL section_vals_val_get(tmp_section, "RHO_ZERO", &
                                         r_val=dft_control%sccs_control%rho_zero)
               CALL cite_reference(fattebert2002)
            CASE DEFAULT
               cpabort("Invalid SCCS model specified. Please, check your input!")
            END SELECT
            CALL cite_reference(yin2017)
         END IF
      END IF
 
      ! Read the planar counter charge section
      tmp_section => section_vals_get_subs_vals(dft_section, "PLANAR_COUNTER_CHARGE")
      CALL section_vals_get(tmp_section, explicit=is_present)
      IF (is_present) THEN
         ! Check section parameter if planar counter charge is activated
         CALL section_vals_val_get(tmp_section, "_SECTION_PARAMETERS_", &
                                   l_val=dft_control%do_pcc)
         IF (dft_control%do_pcc) THEN
            ALLOCATE (dft_control%pcc_control)
            CALL section_vals_val_get(tmp_section, "DIST_EDGE", &
                                      r_val=dft_control%pcc_control%dist_edge)
            CALL section_vals_val_get(tmp_section, "GAU_C", &
                                      r_val=dft_control%pcc_control%gau_c)
            CALL section_vals_val_get(tmp_section, "PARALLEL_PLANE", &
                                      i_val=dft_control%pcc_control%surf_normal)
         END IF
      END IF
 
      ! Read the planar averaged Hartree potential section
      tmp_section => section_vals_get_subs_vals(dft_section, "PLANAR_AVERAGED_V_HARTREE")
      CALL section_vals_get(tmp_section, explicit=is_present)
      IF (is_present) THEN
         ! Check section parameter if the planar-averaged potential is activated
         CALL section_vals_val_get(tmp_section, "_SECTION_PARAMETERS_", &
                                   l_val=dft_control%do_paep)
         IF (dft_control%do_paep) THEN
            ALLOCATE (dft_control%paep_control)
            CALL section_vals_val_get(tmp_section, "PARALLEL_PLANE", &
                                      i_val=dft_control%paep_control%surf_normal)
         END IF
      END IF
 
      ! ZMP added input sections
      ! Read the external density input section
      tmp_section => section_vals_get_subs_vals(dft_section, "EXTERNAL_DENSITY")
      CALL section_vals_get(tmp_section, explicit=dft_control%apply_external_density)
 
      ! Read the external vxc input section
      tmp_section => section_vals_get_subs_vals(dft_section, "EXTERNAL_VXC")
      CALL section_vals_get(tmp_section, explicit=dft_control%apply_external_vxc)
 
      tmp_section => section_vals_get_subs_vals(dft_section, "LOCALIZE")
      CALL section_vals_val_get(tmp_section, "EACH", i_val=dft_control%localize_each)
 
      ! SMEAGOL interface
      tmp_section => section_vals_get_subs_vals(dft_section, "SMEAGOL")
      CALL read_smeagol_control(dft_control%smeagol_control, tmp_section)
 
   END SUBROUTINE read_dft_control
 
! **************************************************************************************************
!> \brief Reads the input and stores in the rixs_control_type
!> \param rixs_control ...
!> \param rixs_section ...
!> \param qs_control ...
! **************************************************************************************************
   SUBROUTINE read_rixs_control(rixs_control, rixs_section, qs_control)
      TYPE(rixs_control_type), POINTER                   :: rixs_control
      TYPE(section_vals_type), POINTER                   :: rixs_section
      TYPE(qs_control_type), POINTER                     :: qs_control
 
      TYPE(section_vals_type), POINTER                   :: td_section, xas_section
 
      CALL section_vals_val_get(rixs_section, "_SECTION_PARAMETERS_", l_val=rixs_control%enabled)
 
      CALL section_vals_val_get(rixs_section, "CORE_STATES", i_val=rixs_control%core_states)
      CALL section_vals_val_get(rixs_section, "VALENCE_STATES", i_val=rixs_control%valence_states)
 
      td_section => section_vals_get_subs_vals(rixs_section, "TDDFPT")
      CALL read_tddfpt2_control(rixs_control%tddfpt2_control, td_section, qs_control)
 
      xas_section => section_vals_get_subs_vals(rixs_section, "XAS_TDP")
      CALL read_xas_tdp_control(rixs_control%xas_tdp_control, xas_section)
 
   END SUBROUTINE read_rixs_control
 
! **************************************************************************************************
!> \brief ...
!> \param qs_control ...
!> \param dft_section ...
! **************************************************************************************************
   SUBROUTINE read_mgrid_section(qs_control, dft_section)
 
      TYPE(qs_control_type), INTENT(INOUT)               :: qs_control
      TYPE(section_vals_type), POINTER                   :: dft_section
 
      CHARACTER(len=*), PARAMETER :: routinen = 'read_mgrid_section'
 
      INTEGER                                            :: handle, igrid_level, ngrid_level
      LOGICAL                                            :: explicit, multigrid_set
      REAL(dp)                                           :: cutoff
      REAL(dp), DIMENSION(:), POINTER                    :: cutofflist
      TYPE(section_vals_type), POINTER                   :: mgrid_section
 
      CALL timeset(routinen, handle)
 
      NULLIFY (mgrid_section, cutofflist)
      mgrid_section => section_vals_get_subs_vals(dft_section, "MGRID")
 
      CALL section_vals_val_get(mgrid_section, "NGRIDS", i_val=ngrid_level)
      CALL section_vals_val_get(mgrid_section, "MULTIGRID_SET", l_val=multigrid_set)
      CALL section_vals_val_get(mgrid_section, "CUTOFF", r_val=cutoff)
      CALL section_vals_val_get(mgrid_section, "PROGRESSION_FACTOR", r_val=qs_control%progression_factor)
      CALL section_vals_val_get(mgrid_section, "COMMENSURATE", l_val=qs_control%commensurate_mgrids)
      CALL section_vals_val_get(mgrid_section, "REALSPACE", l_val=qs_control%realspace_mgrids)
      CALL section_vals_val_get(mgrid_section, "REL_CUTOFF", r_val=qs_control%relative_cutoff)
      CALL section_vals_val_get(mgrid_section, "SKIP_LOAD_BALANCE_DISTRIBUTED", &
                                l_val=qs_control%skip_load_balance_distributed)
 
      ! For SE and DFTB possibly override with new defaults
      IF (qs_control%semi_empirical .OR. qs_control%dftb .OR. qs_control%xtb) THEN
         ngrid_level = 1
         multigrid_set = .false.
         ! Override default cutoff value unless user specified an explicit argument..
         CALL section_vals_val_get(mgrid_section, "CUTOFF", explicit=explicit, r_val=cutoff)
         IF (.NOT. explicit) cutoff = 1.0_dp
      END IF
 
      ALLOCATE (qs_control%e_cutoff(ngrid_level))
      qs_control%cutoff = cutoff
 
      IF (multigrid_set) THEN
         ! Read the values from input
         IF (qs_control%commensurate_mgrids) THEN
            cpabort("Do not specify cutoffs for the commensurate grids (NYI)")
         END IF
 
         CALL section_vals_val_get(mgrid_section, "MULTIGRID_CUTOFF", r_vals=cutofflist)
         IF (ASSOCIATED(cutofflist)) THEN
            IF (SIZE(cutofflist, 1) /= ngrid_level) THEN
               cpabort("Number of multi-grids requested and number of cutoff values do not match")
            END IF
            DO igrid_level = 1, ngrid_level
               qs_control%e_cutoff(igrid_level) = cutofflist(igrid_level)
            END DO
         END IF
         ! set cutoff to smallest value in multgrid available with >= cutoff
         DO igrid_level = ngrid_level, 1, -1
            IF (qs_control%cutoff <= qs_control%e_cutoff(igrid_level)) THEN
               qs_control%cutoff = qs_control%e_cutoff(igrid_level)
               EXIT
            END IF
            ! set largest grid value to cutoff
            IF (igrid_level == 1) THEN
               qs_control%cutoff = qs_control%e_cutoff(1)
            END IF
         END DO
      ELSE
         IF (qs_control%commensurate_mgrids) qs_control%progression_factor = 4.0_dp
         qs_control%e_cutoff(1) = qs_control%cutoff
         DO igrid_level = 2, ngrid_level
            qs_control%e_cutoff(igrid_level) = qs_control%e_cutoff(igrid_level - 1)/ &
                                               qs_control%progression_factor
         END DO
      END IF
      ! check that multigrids are ordered
      DO igrid_level = 2, ngrid_level
         IF (qs_control%e_cutoff(igrid_level) > qs_control%e_cutoff(igrid_level - 1)) THEN
            cpabort("The cutoff values for the multi-grids are not ordered from large to small")
         ELSE IF (qs_control%e_cutoff(igrid_level) == qs_control%e_cutoff(igrid_level - 1)) THEN
            cpabort("The same cutoff value was specified for two multi-grids")
         END IF
      END DO
      CALL timestop(handle)
   END SUBROUTINE read_mgrid_section
 
! **************************************************************************************************
!> \brief ...
!> \param qs_control ...
!> \param qs_section ...
!> \param cell optional cell used to transform Cartesian input vectors
! **************************************************************************************************
   SUBROUTINE read_qs_section(qs_control, qs_section, cell)
 
      TYPE(qs_control_type), INTENT(INOUT)               :: qs_control
      TYPE(section_vals_type), POINTER                   :: qs_section
      TYPE(cell_type), OPTIONAL, POINTER                 :: cell
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'read_qs_section'
 
      CHARACTER(LEN=default_string_length)               :: cval
      CHARACTER(LEN=default_string_length), &
         DIMENSION(:), POINTER                           :: clist
      INTEGER                                            :: handle, itmp, j, jj, k, n_rep, n_var, &
                                                            ngauss, ngp, nrep
      INTEGER, DIMENSION(:), POINTER                     :: tmplist
      LOGICAL :: dftb_scc_mixer_explicit, dftb_tblite_mixer_explicit, explicit, &
         tblite_reference_cli, tblite_reference_cli_section, tblite_section_active, was_present, &
         xtb_scc_mixer_explicit, xtb_tblite_mixer_explicit
      REAL(dp)                                           :: tmp, tmpsqrt, value
      REAL(dp), POINTER                                  :: scal(:)
      TYPE(section_vals_type), POINTER :: cdft_control_section, ddapc_restraint_section, &
         dftb_parameter, dftb_section, dftb_tblite_mixer, eeq_section, genpot_section, &
         lri_optbas_section, mull_section, nonbonded_section, s2_restraint_section, se_section, &
         xtb_parameter, xtb_section, xtb_tblite, xtb_tblite_mixer, xtb_tblite_ref_cli
 
      CALL timeset(routinen, handle)
 
      was_present = .false.
      NULLIFY (mull_section, ddapc_restraint_section, s2_restraint_section, &
               se_section, dftb_section, xtb_section, dftb_parameter, xtb_parameter, lri_optbas_section, &
               cdft_control_section, genpot_section, eeq_section, dftb_tblite_mixer, &
               xtb_tblite_mixer, xtb_tblite_ref_cli)
 
      mull_section => section_vals_get_subs_vals(qs_section, "MULLIKEN_RESTRAINT")
      ddapc_restraint_section => section_vals_get_subs_vals(qs_section, "DDAPC_RESTRAINT")
      s2_restraint_section => section_vals_get_subs_vals(qs_section, "S2_RESTRAINT")
      se_section => section_vals_get_subs_vals(qs_section, "SE")
      dftb_section => section_vals_get_subs_vals(qs_section, "DFTB")
      xtb_section => section_vals_get_subs_vals(qs_section, "xTB")
      dftb_parameter => section_vals_get_subs_vals(dftb_section, "PARAMETER")
      dftb_tblite_mixer => section_vals_get_subs_vals(dftb_section, "TBLITE_MIXER")
      xtb_parameter => section_vals_get_subs_vals(xtb_section, "PARAMETER")
      eeq_section => section_vals_get_subs_vals(xtb_section, "EEQ")
      lri_optbas_section => section_vals_get_subs_vals(qs_section, "OPTIMIZE_LRI_BASIS")
      cdft_control_section => section_vals_get_subs_vals(qs_section, "CDFT")
      nonbonded_section => section_vals_get_subs_vals(xtb_section, "NONBONDED")
      genpot_section => section_vals_get_subs_vals(nonbonded_section, "GENPOT")
      xtb_tblite_mixer => section_vals_get_subs_vals(xtb_section, "TBLITE_MIXER")
      xtb_tblite => section_vals_get_subs_vals(xtb_section, "TBLITE")
      xtb_tblite_ref_cli => section_vals_get_subs_vals(xtb_tblite, "REFERENCE_CLI")
 
      ! Setup all defaults values and overwrite input parameters
      ! EPS_DEFAULT should set the target accuracy in the total energy (~per electron) or a closely related value
      CALL section_vals_val_get(qs_section, "EPS_DEFAULT", r_val=value)
      tmpsqrt = sqrt(value) ! a trick to work around a NAG 5.1 optimizer bug
 
      ! random choice ?
      qs_control%eps_core_charge = value/100.0_dp
      ! correct if all Gaussians would have the same radius (overlap will be smaller than eps_pgf_orb**2).
      ! Can be significantly in error if not... requires fully new screening/pairlist procedures
      qs_control%eps_pgf_orb = tmpsqrt
      qs_control%eps_kg_orb = qs_control%eps_pgf_orb
      ! consistent since also a kind of overlap
      qs_control%eps_ppnl = qs_control%eps_pgf_orb/100.0_dp
      ! accuracy is basically set by the overlap, this sets an empirical shift
      qs_control%eps_ppl = 1.0e-2_dp
      !
      qs_control%gapw_control%eps_cpc = value
      ! expexted error in the density
      qs_control%eps_rho_gspace = value
      qs_control%eps_rho_rspace = value
      ! error in the gradient, can be the sqrt of the error in the energy, ignored if map_consistent
      qs_control%eps_gvg_rspace = tmpsqrt
      !
      CALL section_vals_val_get(qs_section, "EPS_CORE_CHARGE", n_rep_val=n_rep)
      IF (n_rep /= 0) THEN
         CALL section_vals_val_get(qs_section, "EPS_CORE_CHARGE", r_val=qs_control%eps_core_charge)
      END IF
      CALL section_vals_val_get(qs_section, "EPS_GVG_RSPACE", n_rep_val=n_rep)
      IF (n_rep /= 0) THEN
         CALL section_vals_val_get(qs_section, "EPS_GVG_RSPACE", r_val=qs_control%eps_gvg_rspace)
      END IF
      CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", n_rep_val=n_rep)
      IF (n_rep /= 0) THEN
         CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", r_val=qs_control%eps_pgf_orb)
      END IF
      CALL section_vals_val_get(qs_section, "EPS_KG_ORB", n_rep_val=n_rep)
      IF (n_rep /= 0) THEN
         CALL section_vals_val_get(qs_section, "EPS_KG_ORB", r_val=tmp)
         qs_control%eps_kg_orb = sqrt(tmp)
      END IF
      CALL section_vals_val_get(qs_section, "EPS_PPL", n_rep_val=n_rep)
      IF (n_rep /= 0) THEN
         CALL section_vals_val_get(qs_section, "EPS_PPL", r_val=qs_control%eps_ppl)
      END IF
      CALL section_vals_val_get(qs_section, "EPS_PPNL", n_rep_val=n_rep)
      IF (n_rep /= 0) THEN
         CALL section_vals_val_get(qs_section, "EPS_PPNL", r_val=qs_control%eps_ppnl)
      END IF
      CALL section_vals_val_get(qs_section, "EPS_RHO", n_rep_val=n_rep)
      IF (n_rep /= 0) THEN
         CALL section_vals_val_get(qs_section, "EPS_RHO", r_val=qs_control%eps_rho_gspace)
         qs_control%eps_rho_rspace = qs_control%eps_rho_gspace
      END IF
      CALL section_vals_val_get(qs_section, "EPS_RHO_RSPACE", n_rep_val=n_rep)
      IF (n_rep /= 0) THEN
         CALL section_vals_val_get(qs_section, "EPS_RHO_RSPACE", r_val=qs_control%eps_rho_rspace)
      END IF
      CALL section_vals_val_get(qs_section, "EPS_RHO_GSPACE", n_rep_val=n_rep)
      IF (n_rep /= 0) THEN
         CALL section_vals_val_get(qs_section, "EPS_RHO_GSPACE", r_val=qs_control%eps_rho_gspace)
      END IF
      CALL section_vals_val_get(qs_section, "EPS_FILTER_MATRIX", n_rep_val=n_rep)
      IF (n_rep /= 0) THEN
         CALL section_vals_val_get(qs_section, "EPS_FILTER_MATRIX", r_val=qs_control%eps_filter_matrix)
      END IF
      CALL section_vals_val_get(qs_section, "EPS_CPC", n_rep_val=n_rep)
      IF (n_rep /= 0) THEN
         CALL section_vals_val_get(qs_section, "EPS_CPC", r_val=qs_control%gapw_control%eps_cpc)
      END IF
 
      CALL section_vals_val_get(qs_section, "EPSFIT", r_val=qs_control%gapw_control%eps_fit)
      CALL section_vals_val_get(qs_section, "EPSISO", r_val=qs_control%gapw_control%eps_iso)
      CALL section_vals_val_get(qs_section, "EPSSVD", r_val=qs_control%gapw_control%eps_svd)
      CALL section_vals_val_get(qs_section, "EPSRHO0", r_val=qs_control%gapw_control%eps_Vrho0)
      CALL section_vals_val_get(qs_section, "ALPHA0_HARD", r_val=qs_control%gapw_control%alpha0_hard)
      qs_control%gapw_control%alpha0_hard_from_input = .false.
      IF (qs_control%gapw_control%alpha0_hard /= 0.0_dp) qs_control%gapw_control%alpha0_hard_from_input = .true.
      CALL section_vals_val_get(qs_section, "FORCE_PAW", l_val=qs_control%gapw_control%force_paw)
      CALL section_vals_val_get(qs_section, "MAX_RAD_LOCAL", r_val=qs_control%gapw_control%max_rad_local)
 
      CALL section_vals_val_get(qs_section, "MIN_PAIR_LIST_RADIUS", r_val=qs_control%pairlist_radius)
 
      CALL section_vals_val_get(qs_section, "LS_SCF", l_val=qs_control%do_ls_scf)
      CALL section_vals_val_get(qs_section, "ALMO_SCF", l_val=qs_control%do_almo_scf)
      CALL section_vals_val_get(qs_section, "KG_METHOD", l_val=qs_control%do_kg)
 
      ! Logicals
      CALL section_vals_val_get(qs_section, "REF_EMBED_SUBSYS", l_val=qs_control%ref_embed_subsys)
      CALL section_vals_val_get(qs_section, "CLUSTER_EMBED_SUBSYS", l_val=qs_control%cluster_embed_subsys)
      CALL section_vals_val_get(qs_section, "HIGH_LEVEL_EMBED_SUBSYS", l_val=qs_control%high_level_embed_subsys)
      CALL section_vals_val_get(qs_section, "DFET_EMBEDDED", l_val=qs_control%dfet_embedded)
      CALL section_vals_val_get(qs_section, "DMFET_EMBEDDED", l_val=qs_control%dmfet_embedded)
 
      ! Integers gapw
      CALL section_vals_val_get(qs_section, "LMAXN1", i_val=qs_control%gapw_control%lmax_sphere)
      CALL section_vals_val_get(qs_section, "LMAXN0", i_val=qs_control%gapw_control%lmax_rho0)
      CALL section_vals_val_get(qs_section, "LADDN0", i_val=qs_control%gapw_control%ladd_rho0)
      CALL section_vals_val_get(qs_section, "QUADRATURE", i_val=qs_control%gapw_control%quadrature)
      ! GAPW 1c basis
      CALL section_vals_val_get(qs_section, "GAPW_1C_BASIS", i_val=qs_control%gapw_control%basis_1c)
      IF (qs_control%gapw_control%basis_1c /= gapw_1c_orb) THEN
         qs_control%gapw_control%eps_svd = max(qs_control%gapw_control%eps_svd, 1.e-12_dp)
      END IF
      ! GAPW accurate integration
      CALL section_vals_val_get(qs_section, "GAPW_ACCURATE_XCINT", l_val=qs_control%gapw_control%accurate_xcint)
      CALL section_vals_val_get(qs_section, "ALPHA_WEIGHTS", r_val=qs_control%gapw_control%aweights)
 
      ! Integers grids
      CALL section_vals_val_get(qs_section, "PW_GRID", i_val=itmp)
      SELECT CASE (itmp)
      CASE (do_pwgrid_spherical)
         qs_control%pw_grid_opt%spherical = .true.
         qs_control%pw_grid_opt%fullspace = .false.
      CASE (do_pwgrid_ns_fullspace)
         qs_control%pw_grid_opt%spherical = .false.
         qs_control%pw_grid_opt%fullspace = .true.
      CASE (do_pwgrid_ns_halfspace)
         qs_control%pw_grid_opt%spherical = .false.
         qs_control%pw_grid_opt%fullspace = .false.
      END SELECT
 
      !   Method for PPL calculation
      CALL section_vals_val_get(qs_section, "CORE_PPL", i_val=itmp)
      qs_control%do_ppl_method = itmp
 
      CALL section_vals_val_get(qs_section, "PW_GRID_LAYOUT", i_vals=tmplist)
      qs_control%pw_grid_opt%distribution_layout = tmplist
      CALL section_vals_val_get(qs_section, "PW_GRID_BLOCKED", i_val=qs_control%pw_grid_opt%blocked)
 
      !Integers extrapolation
      CALL section_vals_val_get(qs_section, "EXTRAPOLATION", i_val=qs_control%wf_interpolation_method_nr)
      CALL section_vals_val_get(qs_section, "EXTRAPOLATION_ORDER", i_val=qs_control%wf_extrapolation_order)
 
      !Method
      CALL section_vals_val_get(qs_section, "METHOD", i_val=qs_control%method_id)
      qs_control%gapw = .false.
      qs_control%gapw_xc = .false.
      qs_control%gpw = .false.
      qs_control%pao = .false.
      qs_control%dftb = .false.
      qs_control%xtb = .false.
      qs_control%semi_empirical = .false.
      qs_control%ofgpw = .false.
      qs_control%lrigpw = .false.
      qs_control%rigpw = .false.
      SELECT CASE (qs_control%method_id)
      CASE (do_method_gapw)
         CALL cite_reference(lippert1999)
         CALL cite_reference(krack2000)
         qs_control%gapw = .true.
      CASE (do_method_gapw_xc)
         qs_control%gapw_xc = .true.
      CASE (do_method_gpw)
         CALL cite_reference(lippert1997)
         CALL cite_reference(vandevondele2005a)
         qs_control%gpw = .true.
      CASE (do_method_ofgpw)
         qs_control%ofgpw = .true.
      CASE (do_method_lrigpw)
         qs_control%lrigpw = .true.
      CASE (do_method_rigpw)
         qs_control%rigpw = .true.
      CASE (do_method_dftb)
         qs_control%dftb = .true.
         CALL cite_reference(porezag1995)
         CALL cite_reference(seifert1996)
      CASE (do_method_xtb)
         qs_control%xtb = .true.
         CALL cite_reference(grimme2017)
         CALL cite_reference(pracht2019)
      CASE (do_method_mndo)
         CALL cite_reference(dewar1977)
         qs_control%semi_empirical = .true.
      CASE (do_method_am1)
         CALL cite_reference(dewar1985)
         qs_control%semi_empirical = .true.
      CASE (do_method_pm3)
         CALL cite_reference(stewart1989)
         qs_control%semi_empirical = .true.
      CASE (do_method_pnnl)
         CALL cite_reference(schenter2008)
         qs_control%semi_empirical = .true.
      CASE (do_method_pm6)
         CALL cite_reference(stewart2007)
         qs_control%semi_empirical = .true.
      CASE (do_method_pm6fm)
         CALL cite_reference(vanvoorhis2015)
         qs_control%semi_empirical = .true.
      CASE (do_method_pdg)
         CALL cite_reference(repasky2002)
         qs_control%semi_empirical = .true.
      CASE (do_method_rm1)
         CALL cite_reference(rocha2006)
         qs_control%semi_empirical = .true.
      CASE (do_method_mndod)
         CALL cite_reference(dewar1977)
         CALL cite_reference(thiel1992)
         qs_control%semi_empirical = .true.
      END SELECT
 
      CALL section_vals_get(mull_section, explicit=qs_control%mulliken_restraint)
 
      IF (qs_control%mulliken_restraint) THEN
         CALL section_vals_val_get(mull_section, "STRENGTH", r_val=qs_control%mulliken_restraint_control%strength)
         CALL section_vals_val_get(mull_section, "TARGET", r_val=qs_control%mulliken_restraint_control%target)
         CALL section_vals_val_get(mull_section, "ATOMS", n_rep_val=n_rep)
         jj = 0
         DO k = 1, n_rep
            CALL section_vals_val_get(mull_section, "ATOMS", i_rep_val=k, i_vals=tmplist)
            jj = jj + SIZE(tmplist)
         END DO
         qs_control%mulliken_restraint_control%natoms = jj
         IF (qs_control%mulliken_restraint_control%natoms < 1) &
            cpabort("Need at least 1 atom to use mulliken constraints")
         ALLOCATE (qs_control%mulliken_restraint_control%atoms(qs_control%mulliken_restraint_control%natoms))
         jj = 0
         DO k = 1, n_rep
            CALL section_vals_val_get(mull_section, "ATOMS", i_rep_val=k, i_vals=tmplist)
            DO j = 1, SIZE(tmplist)
               jj = jj + 1
               qs_control%mulliken_restraint_control%atoms(jj) = tmplist(j)
            END DO
         END DO
      END IF
      CALL section_vals_get(ddapc_restraint_section, n_repetition=nrep, explicit=qs_control%ddapc_restraint)
      IF (qs_control%ddapc_restraint) THEN
         ALLOCATE (qs_control%ddapc_restraint_control(nrep))
         CALL read_ddapc_section(qs_control, qs_section=qs_section)
         qs_control%ddapc_restraint_is_spin = .false.
         qs_control%ddapc_explicit_potential = .false.
      END IF
 
      CALL section_vals_get(s2_restraint_section, explicit=qs_control%s2_restraint)
      IF (qs_control%s2_restraint) THEN
         CALL section_vals_val_get(s2_restraint_section, "STRENGTH", &
                                   r_val=qs_control%s2_restraint_control%strength)
         CALL section_vals_val_get(s2_restraint_section, "TARGET", &
                                   r_val=qs_control%s2_restraint_control%target)
         CALL section_vals_val_get(s2_restraint_section, "FUNCTIONAL_FORM", &
                                   i_val=qs_control%s2_restraint_control%functional_form)
      END IF
 
      CALL section_vals_get(cdft_control_section, explicit=qs_control%cdft)
      IF (qs_control%cdft) THEN
         CALL read_cdft_control_section(qs_control, cdft_control_section)
      END IF
 
      ! Semi-empirical code
      IF (qs_control%semi_empirical) THEN
         CALL section_vals_val_get(se_section, "ORTHOGONAL_BASIS", &
                                   l_val=qs_control%se_control%orthogonal_basis)
         CALL section_vals_val_get(se_section, "DELTA", &
                                   r_val=qs_control%se_control%delta)
         CALL section_vals_val_get(se_section, "ANALYTICAL_GRADIENTS", &
                                   l_val=qs_control%se_control%analytical_gradients)
         CALL section_vals_val_get(se_section, "FORCE_KDSO-D_EXCHANGE", &
                                   l_val=qs_control%se_control%force_kdsod_EX)
         ! Integral Screening
         CALL section_vals_val_get(se_section, "INTEGRAL_SCREENING", &
                                   i_val=qs_control%se_control%integral_screening)
         IF (qs_control%method_id == do_method_pnnl) THEN
            IF (qs_control%se_control%integral_screening /= do_se_is_slater) &
               CALL cp_warn(__location__, &
                            "PNNL semi-empirical parameterization supports only the Slater type "// &
                            "integral scheme. Revert to Slater and continue the calculation.")
            qs_control%se_control%integral_screening = do_se_is_slater
         END IF
         ! Global Arrays variable
         CALL section_vals_val_get(se_section, "GA%NCELLS", &
                                   i_val=qs_control%se_control%ga_ncells)
         ! Long-Range correction
         CALL section_vals_val_get(se_section, "LR_CORRECTION%CUTOFF", &
                                   r_val=qs_control%se_control%cutoff_lrc)
         qs_control%se_control%taper_lrc = qs_control%se_control%cutoff_lrc
         CALL section_vals_val_get(se_section, "LR_CORRECTION%RC_TAPER", &
                                   explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(se_section, "LR_CORRECTION%RC_TAPER", &
                                      r_val=qs_control%se_control%taper_lrc)
         END IF
         CALL section_vals_val_get(se_section, "LR_CORRECTION%RC_RANGE", &
                                   r_val=qs_control%se_control%range_lrc)
         ! Coulomb
         CALL section_vals_val_get(se_section, "COULOMB%CUTOFF", &
                                   r_val=qs_control%se_control%cutoff_cou)
         qs_control%se_control%taper_cou = qs_control%se_control%cutoff_cou
         CALL section_vals_val_get(se_section, "COULOMB%RC_TAPER", &
                                   explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(se_section, "COULOMB%RC_TAPER", &
                                      r_val=qs_control%se_control%taper_cou)
         END IF
         CALL section_vals_val_get(se_section, "COULOMB%RC_RANGE", &
                                   r_val=qs_control%se_control%range_cou)
         ! Exchange
         CALL section_vals_val_get(se_section, "EXCHANGE%CUTOFF", &
                                   r_val=qs_control%se_control%cutoff_exc)
         qs_control%se_control%taper_exc = qs_control%se_control%cutoff_exc
         CALL section_vals_val_get(se_section, "EXCHANGE%RC_TAPER", &
                                   explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(se_section, "EXCHANGE%RC_TAPER", &
                                      r_val=qs_control%se_control%taper_exc)
         END IF
         CALL section_vals_val_get(se_section, "EXCHANGE%RC_RANGE", &
                                   r_val=qs_control%se_control%range_exc)
         ! Screening (only if the integral scheme is of dumped type)
         IF (qs_control%se_control%integral_screening == do_se_is_kdso_d) THEN
            CALL section_vals_val_get(se_section, "SCREENING%RC_TAPER", &
                                      r_val=qs_control%se_control%taper_scr)
            CALL section_vals_val_get(se_section, "SCREENING%RC_RANGE", &
                                      r_val=qs_control%se_control%range_scr)
         END IF
         ! Periodic Type Calculation
         CALL section_vals_val_get(se_section, "PERIODIC", &
                                   i_val=qs_control%se_control%periodic_type)
         SELECT CASE (qs_control%se_control%periodic_type)
         CASE (do_se_lr_none)
            qs_control%se_control%do_ewald = .false.
            qs_control%se_control%do_ewald_r3 = .false.
            qs_control%se_control%do_ewald_gks = .false.
         CASE (do_se_lr_ewald)
            qs_control%se_control%do_ewald = .true.
            qs_control%se_control%do_ewald_r3 = .false.
            qs_control%se_control%do_ewald_gks = .false.
         CASE (do_se_lr_ewald_gks)
            qs_control%se_control%do_ewald = .false.
            qs_control%se_control%do_ewald_r3 = .false.
            qs_control%se_control%do_ewald_gks = .true.
            IF (qs_control%method_id /= do_method_pnnl) &
               CALL cp_abort(__location__, &
                             "A periodic semi-empirical calculation was requested with a long-range  "// &
                             "summation on the single integral evaluation. This scheme is supported  "// &
                             "only by the PNNL parameterization.")
         CASE (do_se_lr_ewald_r3)
            qs_control%se_control%do_ewald = .true.
            qs_control%se_control%do_ewald_r3 = .true.
            qs_control%se_control%do_ewald_gks = .false.
            IF (qs_control%se_control%integral_screening /= do_se_is_kdso) &
               CALL cp_abort(__location__, &
                             "A periodic semi-empirical calculation was requested with a long-range  "// &
                             "summation for the slowly convergent part 1/R^3, which is not congruent "// &
                             "with the integral screening chosen. The only integral screening supported "// &
                             "by this periodic type calculation is the standard Klopman-Dewar-Sabelli-Ohno.")
         END SELECT
 
         ! dispersion pair potentials
         CALL section_vals_val_get(se_section, "DISPERSION", &
                                   l_val=qs_control%se_control%dispersion)
         CALL section_vals_val_get(se_section, "DISPERSION_RADIUS", &
                                   r_val=qs_control%se_control%rcdisp)
         CALL section_vals_val_get(se_section, "COORDINATION_CUTOFF", &
                                   r_val=qs_control%se_control%epscn)
         CALL section_vals_val_get(se_section, "D3_SCALING", r_vals=scal)
         qs_control%se_control%sd3(1) = scal(1)
         qs_control%se_control%sd3(2) = scal(2)
         qs_control%se_control%sd3(3) = scal(3)
         CALL section_vals_val_get(se_section, "DISPERSION_PARAMETER_FILE", &
                                   c_val=qs_control%se_control%dispersion_parameter_file)
 
         ! Stop the execution for non-implemented features
         IF (qs_control%se_control%periodic_type == do_se_lr_ewald_r3) THEN
            cpabort("EWALD_R3 not implemented yet!")
         END IF
 
         IF (qs_control%method_id == do_method_mndo .OR. &
             qs_control%method_id == do_method_am1 .OR. &
             qs_control%method_id == do_method_mndod .OR. &
             qs_control%method_id == do_method_pdg .OR. &
             qs_control%method_id == do_method_pm3 .OR. &
             qs_control%method_id == do_method_pm6 .OR. &
             qs_control%method_id == do_method_pm6fm .OR. &
             qs_control%method_id == do_method_pnnl .OR. &
             qs_control%method_id == do_method_rm1) THEN
            qs_control%se_control%orthogonal_basis = .true.
         END IF
      END IF
 
      ! DFTB code
      IF (qs_control%dftb) THEN
         CALL section_vals_val_get(dftb_section, "ORTHOGONAL_BASIS", &
                                   l_val=qs_control%dftb_control%orthogonal_basis)
         CALL section_vals_val_get(dftb_section, "SELF_CONSISTENT", &
                                   l_val=qs_control%dftb_control%self_consistent)
         CALL section_vals_val_get(dftb_section, "DISPERSION", &
                                   l_val=qs_control%dftb_control%dispersion)
         CALL section_vals_val_get(dftb_section, "DIAGONAL_DFTB3", &
                                   l_val=qs_control%dftb_control%dftb3_diagonal)
         CALL section_vals_val_get(dftb_section, "HB_SR_GAMMA", &
                                   l_val=qs_control%dftb_control%hb_sr_damp)
         CALL section_vals_val_get(dftb_section, "SCC_MIXER", &
                                   explicit=dftb_scc_mixer_explicit)
         CALL section_vals_val_get(dftb_section, "SCC_MIXER", &
                                   i_val=qs_control%dftb_control%tblite_scc_mixer)
         CALL section_vals_get(dftb_tblite_mixer, explicit=dftb_tblite_mixer_explicit)
         CALL read_tblite_mixer_section(dftb_tblite_mixer, &
                                        qs_control%dftb_control%tblite_mixer_iterations, &
                                        qs_control%dftb_control%tblite_mixer_memory, &
                                        qs_control%dftb_control%tblite_mixer_solver, &
                                        qs_control%dftb_control%tblite_mixer_damping, &
                                        qs_control%dftb_control%tblite_mixer_omega0, &
                                        qs_control%dftb_control%tblite_mixer_min_weight, &
                                        qs_control%dftb_control%tblite_mixer_max_weight, &
                                        qs_control%dftb_control%tblite_mixer_weight_factor, &
                                        "DFTB/TBLITE_MIXER")
         IF (qs_control%do_ls_scf) THEN
            IF (dftb_scc_mixer_explicit .AND. &
                qs_control%dftb_control%tblite_scc_mixer /= tblite_scc_mixer_none) &
               CALL cp_warn(__location__, &
                            "DFTB/SCC_MIXER is reset to NONE with QS/LS_SCF; LS_SCF optimizes "// &
                            "the density matrix directly.")
            IF (dftb_tblite_mixer_explicit) &
               CALL cp_warn(__location__, &
                            "DFTB/TBLITE_MIXER settings are ignored with QS/LS_SCF; LS_SCF controls "// &
                            "the density-matrix optimization.")
            qs_control%dftb_control%tblite_scc_mixer = tblite_scc_mixer_none
         END IF
         IF (qs_control%dftb_control%tblite_mixer_damping <= 0.0_dp) &
            cpabort("DFTB/TBLITE_MIXER/DAMPING must be positive")
         CALL section_vals_val_get(dftb_section, "EPS_DISP", &
                                   r_val=qs_control%dftb_control%eps_disp)
         CALL section_vals_val_get(dftb_section, "DO_EWALD", explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(dftb_section, "DO_EWALD", &
                                      l_val=qs_control%dftb_control%do_ewald)
         ELSE
            qs_control%dftb_control%do_ewald = (qs_control%periodicity /= 0)
         END IF
         CALL section_vals_val_get(dftb_parameter, "PARAM_FILE_PATH", &
                                   c_val=qs_control%dftb_control%sk_file_path)
         CALL section_vals_val_get(dftb_parameter, "PARAM_FILE_NAME", &
                                   c_val=qs_control%dftb_control%sk_file_list)
         CALL section_vals_val_get(dftb_parameter, "HB_SR_PARAM", &
                                   r_val=qs_control%dftb_control%hb_sr_para)
         CALL section_vals_val_get(dftb_parameter, "SK_FILE", n_rep_val=n_var)
         ALLOCATE (qs_control%dftb_control%sk_pair_list(3, n_var))
         DO k = 1, n_var
            CALL section_vals_val_get(dftb_parameter, "SK_FILE", i_rep_val=k, &
                                      c_vals=clist)
            qs_control%dftb_control%sk_pair_list(1:3, k) = clist(1:3)
         END DO
         ! Dispersion type
         CALL section_vals_val_get(dftb_parameter, "DISPERSION_TYPE", &
                                   i_val=qs_control%dftb_control%dispersion_type)
         CALL section_vals_val_get(dftb_parameter, "UFF_FORCE_FIELD", &
                                   c_val=qs_control%dftb_control%uff_force_field)
         ! D3 Dispersion
         CALL section_vals_val_get(dftb_parameter, "DISPERSION_RADIUS", &
                                   r_val=qs_control%dftb_control%rcdisp)
         CALL section_vals_val_get(dftb_parameter, "COORDINATION_CUTOFF", &
                                   r_val=qs_control%dftb_control%epscn)
         CALL section_vals_val_get(dftb_parameter, "D2_EXP_PRE", &
                                   r_val=qs_control%dftb_control%exp_pre)
         CALL section_vals_val_get(dftb_parameter, "D2_SCALING", &
                                   r_val=qs_control%dftb_control%scaling)
         CALL section_vals_val_get(dftb_parameter, "D3_SCALING", r_vals=scal)
         qs_control%dftb_control%sd3(1) = scal(1)
         qs_control%dftb_control%sd3(2) = scal(2)
         qs_control%dftb_control%sd3(3) = scal(3)
         CALL section_vals_val_get(dftb_parameter, "D3BJ_SCALING", r_vals=scal)
         qs_control%dftb_control%sd3bj(1) = scal(1)
         qs_control%dftb_control%sd3bj(2) = scal(2)
         qs_control%dftb_control%sd3bj(3) = scal(3)
         qs_control%dftb_control%sd3bj(4) = scal(4)
         CALL section_vals_val_get(dftb_parameter, "DISPERSION_PARAMETER_FILE", &
                                   c_val=qs_control%dftb_control%dispersion_parameter_file)
 
         IF (qs_control%dftb_control%dispersion) CALL cite_reference(zhechkov2005)
         IF (qs_control%dftb_control%self_consistent) CALL cite_reference(elstner1998)
         IF (qs_control%dftb_control%hb_sr_damp) CALL cite_reference(hu2007)
      END IF
 
      ! xTB code
      IF (qs_control%xtb) THEN
         CALL section_vals_val_get(xtb_section, "GFN_TYPE", i_val=qs_control%xtb_control%gfn_type)
         CALL section_vals_val_get(xtb_tblite, "_SECTION_PARAMETERS_", l_val=tblite_section_active)
         qs_control%xtb_control%do_tblite = (qs_control%xtb_control%gfn_type == gfn_tblite)
         IF (qs_control%xtb_control%do_tblite) THEN
            IF (.NOT. tblite_section_active) &
               cpabort("XTB/GFN_TYPE TBLITE requires an XTB/TBLITE section")
            ! The CP2K-internal GFN1 defaults are still used to initialize shared xTB fields.
            qs_control%xtb_control%gfn_type = gfn1xtb
         ELSE IF (tblite_section_active) THEN
            cpabort("The XTB/TBLITE section requires XTB/GFN_TYPE TBLITE")
         END IF
         CALL section_vals_val_get(xtb_section, "SCC_MIXER", &
                                   explicit=xtb_scc_mixer_explicit)
         CALL section_vals_val_get(xtb_section, "SCC_MIXER", &
                                   i_val=qs_control%xtb_control%tblite_scc_mixer)
         CALL section_vals_get(xtb_tblite_mixer, explicit=xtb_tblite_mixer_explicit)
         CALL read_tblite_mixer_section(xtb_tblite_mixer, &
                                        qs_control%xtb_control%tblite_mixer_iterations, &
                                        qs_control%xtb_control%tblite_mixer_memory, &
                                        qs_control%xtb_control%tblite_mixer_solver, &
                                        qs_control%xtb_control%tblite_mixer_damping, &
                                        qs_control%xtb_control%tblite_mixer_omega0, &
                                        qs_control%xtb_control%tblite_mixer_min_weight, &
                                        qs_control%xtb_control%tblite_mixer_max_weight, &
                                        qs_control%xtb_control%tblite_mixer_weight_factor, &
                                        "XTB/TBLITE_MIXER")
         IF (xtb_tblite_mixer_explicit) &
            CALL section_vals_val_get(xtb_tblite_mixer, "DAMPING", &
                                      explicit=qs_control%xtb_control%tblite_mixer_damping_explicit)
         IF ((.NOT. qs_control%xtb_control%do_tblite) .AND. &
             qs_control%xtb_control%gfn_type == 0) THEN
            IF (xtb_scc_mixer_explicit .AND. &
                qs_control%xtb_control%tblite_scc_mixer /= tblite_scc_mixer_auto .AND. &
                qs_control%xtb_control%tblite_scc_mixer /= tblite_scc_mixer_none) &
               CALL cp_warn(__location__, &
                            "XTB/SCC_MIXER is reset to NONE for CP2K-internal GFN0-xTB; "// &
                            "GFN0-xTB has no SCC variables to mix.")
            IF (xtb_tblite_mixer_explicit) &
               CALL cp_warn(__location__, &
                            "XTB/TBLITE_MIXER settings are ignored for CP2K-internal GFN0-xTB; "// &
                            "GFN0-xTB has no SCC variables to mix.")
            qs_control%xtb_control%tblite_scc_mixer = tblite_scc_mixer_none
         END IF
         IF (qs_control%do_ls_scf) THEN
            IF (xtb_scc_mixer_explicit .AND. &
                qs_control%xtb_control%tblite_scc_mixer /= tblite_scc_mixer_none) &
               CALL cp_warn(__location__, &
                            "XTB/SCC_MIXER is reset to NONE with QS/LS_SCF; LS_SCF optimizes "// &
                            "the density matrix directly.")
            IF (xtb_tblite_mixer_explicit) &
               CALL cp_warn(__location__, &
                            "XTB/TBLITE_MIXER settings are ignored with QS/LS_SCF; LS_SCF controls "// &
                            "the density-matrix optimization.")
            qs_control%xtb_control%tblite_scc_mixer = tblite_scc_mixer_none
         END IF
         IF (qs_control%xtb_control%tblite_mixer_damping <= 0.0_dp) &
            cpabort("XTB/TBLITE_MIXER/DAMPING must be positive")
         CALL section_vals_val_get(xtb_section, "DO_EWALD", explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(xtb_section, "DO_EWALD", &
                                      l_val=qs_control%xtb_control%do_ewald)
         ELSE
            qs_control%xtb_control%do_ewald = (qs_control%periodicity /= 0)
         END IF
         ! vdW
         CALL section_vals_val_get(xtb_section, "VDW_POTENTIAL", explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(xtb_section, "VDW_POTENTIAL", c_val=cval)
            CALL uppercase(cval)
            SELECT CASE (cval)
            CASE ("NONE")
               qs_control%xtb_control%vdw_type = xtb_vdw_type_none
            CASE ("DFTD3")
               qs_control%xtb_control%vdw_type = xtb_vdw_type_d3
            CASE ("DFTD4")
               qs_control%xtb_control%vdw_type = xtb_vdw_type_d4
            CASE DEFAULT
               cpabort("vdW type")
            END SELECT
         ELSE
            SELECT CASE (qs_control%xtb_control%gfn_type)
            CASE (0)
               qs_control%xtb_control%vdw_type = xtb_vdw_type_d4
            CASE (1)
               qs_control%xtb_control%vdw_type = xtb_vdw_type_d3
            CASE (2)
               qs_control%xtb_control%vdw_type = xtb_vdw_type_d4
               cpabort("gfn2-xtb tbd")
            CASE DEFAULT
               cpabort("GFN type")
            END SELECT
         END IF
         !
         CALL section_vals_val_get(xtb_section, "STO_NG", i_val=ngauss)
         qs_control%xtb_control%sto_ng = ngauss
         CALL section_vals_val_get(xtb_section, "HYDROGEN_STO_NG", i_val=ngauss)
         qs_control%xtb_control%h_sto_ng = ngauss
         CALL section_vals_val_get(xtb_parameter, "PARAM_FILE_PATH", &
                                   c_val=qs_control%xtb_control%parameter_file_path)
         CALL section_vals_val_get(xtb_parameter, "PARAM_FILE_NAME", explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(xtb_parameter, "PARAM_FILE_NAME", &
                                      c_val=qs_control%xtb_control%parameter_file_name)
         ELSE
            SELECT CASE (qs_control%xtb_control%gfn_type)
            CASE (0)
               qs_control%xtb_control%parameter_file_name = "xTB0_parameters"
            CASE (1)
               qs_control%xtb_control%parameter_file_name = "xTB1_parameters"
            CASE (2)
               cpabort("gfn2-xtb tbd")
            CASE DEFAULT
               cpabort("GFN type")
            END SELECT
         END IF
         ! D3 Dispersion
         CALL section_vals_val_get(xtb_parameter, "DISPERSION_RADIUS", &
                                   r_val=qs_control%xtb_control%rcdisp)
         CALL section_vals_val_get(xtb_parameter, "COORDINATION_CUTOFF", &
                                   r_val=qs_control%xtb_control%epscn)
         CALL section_vals_val_get(xtb_parameter, "D3BJ_SCALING", explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(xtb_parameter, "D3BJ_SCALING", r_vals=scal)
            qs_control%xtb_control%s6 = scal(1)
            qs_control%xtb_control%s8 = scal(2)
         ELSE
            SELECT CASE (qs_control%xtb_control%gfn_type)
            CASE (0)
               qs_control%xtb_control%s6 = 1.00_dp
               qs_control%xtb_control%s8 = 2.85_dp
            CASE (1)
               qs_control%xtb_control%s6 = 1.00_dp
               qs_control%xtb_control%s8 = 2.40_dp
            CASE (2)
               cpabort("gfn2-xtb tbd")
            CASE DEFAULT
               cpabort("GFN type")
            END SELECT
         END IF
         CALL section_vals_val_get(xtb_parameter, "D3BJ_PARAM", explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(xtb_parameter, "D3BJ_PARAM", r_vals=scal)
            qs_control%xtb_control%a1 = scal(1)
            qs_control%xtb_control%a2 = scal(2)
         ELSE
            SELECT CASE (qs_control%xtb_control%gfn_type)
            CASE (0)
               qs_control%xtb_control%a1 = 0.80_dp
               qs_control%xtb_control%a2 = 4.60_dp
            CASE (1)
               qs_control%xtb_control%a1 = 0.63_dp
               qs_control%xtb_control%a2 = 5.00_dp
            CASE (2)
               cpabort("gfn2-xtb tbd")
            CASE DEFAULT
               cpabort("GFN type")
            END SELECT
         END IF
         CALL section_vals_val_get(xtb_parameter, "DISPERSION_PARAMETER_FILE", &
                                   c_val=qs_control%xtb_control%dispersion_parameter_file)
         ! global parameters
         CALL section_vals_val_get(xtb_parameter, "HUCKEL_CONSTANTS", explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(xtb_parameter, "HUCKEL_CONSTANTS", r_vals=scal)
            qs_control%xtb_control%ks = scal(1)
            qs_control%xtb_control%kp = scal(2)
            qs_control%xtb_control%kd = scal(3)
            qs_control%xtb_control%ksp = scal(4)
            qs_control%xtb_control%k2sh = scal(5)
            IF (qs_control%xtb_control%gfn_type == 0) THEN
               ! enforce ksp for gfn0
               qs_control%xtb_control%ksp = 0.5_dp*(scal(1) + scal(2))
            END IF
         ELSE
            SELECT CASE (qs_control%xtb_control%gfn_type)
            CASE (0)
               qs_control%xtb_control%ks = 2.00_dp
               qs_control%xtb_control%kp = 2.4868_dp
               qs_control%xtb_control%kd = 2.27_dp
               qs_control%xtb_control%ksp = 2.2434_dp
               qs_control%xtb_control%k2sh = 1.1241_dp
            CASE (1)
               qs_control%xtb_control%ks = 1.85_dp
               qs_control%xtb_control%kp = 2.25_dp
               qs_control%xtb_control%kd = 2.00_dp
               qs_control%xtb_control%ksp = 2.08_dp
               qs_control%xtb_control%k2sh = 2.85_dp
            CASE (2)
               cpabort("gfn2-xtb tbd")
            CASE DEFAULT
               cpabort("GFN type")
            END SELECT
         END IF
         CALL section_vals_val_get(xtb_parameter, "COULOMB_CONSTANTS", explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(xtb_parameter, "COULOMB_CONSTANTS", r_vals=scal)
            qs_control%xtb_control%kg = scal(1)
            qs_control%xtb_control%kf = scal(2)
         ELSE
            SELECT CASE (qs_control%xtb_control%gfn_type)
            CASE (0)
               qs_control%xtb_control%kg = 2.00_dp
               qs_control%xtb_control%kf = 1.50_dp
            CASE (1)
               qs_control%xtb_control%kg = 2.00_dp
               qs_control%xtb_control%kf = 1.50_dp
            CASE (2)
               cpabort("gfn2-xtb tbd")
            CASE DEFAULT
               cpabort("GFN type")
            END SELECT
         END IF
         CALL section_vals_val_get(xtb_parameter, "CN_CONSTANTS", r_vals=scal)
         qs_control%xtb_control%kcns = scal(1)
         qs_control%xtb_control%kcnp = scal(2)
         qs_control%xtb_control%kcnd = scal(3)
         !
         CALL section_vals_val_get(xtb_parameter, "EN_CONSTANTS", explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(xtb_parameter, "EN_CONSTANTS", r_vals=scal)
            SELECT CASE (qs_control%xtb_control%gfn_type)
            CASE (0)
               qs_control%xtb_control%ksen = scal(1)
               qs_control%xtb_control%kpen = scal(2)
               qs_control%xtb_control%kden = scal(3)
            CASE (1)
               qs_control%xtb_control%ken = scal(1)
            CASE (2)
               cpabort("gfn2-xtb tbd")
            CASE DEFAULT
               cpabort("GFN type")
            END SELECT
         ELSE
            SELECT CASE (qs_control%xtb_control%gfn_type)
            CASE (0)
               qs_control%xtb_control%ksen = 0.006_dp
               qs_control%xtb_control%kpen = -0.001_dp
               qs_control%xtb_control%kden = -0.002_dp
            CASE (1)
               qs_control%xtb_control%ken = -0.007_dp
            CASE (2)
               cpabort("gfn2-xtb tbd")
            CASE DEFAULT
               cpabort("GFN type")
            END SELECT
         END IF
         ! ben
         CALL section_vals_val_get(xtb_parameter, "BEN_CONSTANT", r_vals=scal)
         qs_control%xtb_control%ben = scal(1)
         ! enscale (hidden parameter in repulsion
         CALL section_vals_val_get(xtb_parameter, "ENSCALE", explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(xtb_parameter, "ENSCALE", &
                                      r_val=qs_control%xtb_control%enscale)
         ELSE
            SELECT CASE (qs_control%xtb_control%gfn_type)
            CASE (0)
               qs_control%xtb_control%enscale = -0.09_dp
            CASE (1)
               qs_control%xtb_control%enscale = 0._dp
            CASE (2)
               cpabort("gfn2-xtb tbd")
            CASE DEFAULT
               cpabort("GFN type")
            END SELECT
         END IF
         ! XB
         CALL section_vals_val_get(xtb_section, "USE_HALOGEN_CORRECTION", &
                                   l_val=qs_control%xtb_control%xb_interaction)
         CALL section_vals_val_get(xtb_parameter, "HALOGEN_BINDING", r_vals=scal)
         qs_control%xtb_control%kxr = scal(1)
         qs_control%xtb_control%kx2 = scal(2)
         ! NONBONDED interactions
         CALL section_vals_val_get(xtb_section, "DO_NONBONDED", &
                                   l_val=qs_control%xtb_control%do_nonbonded)
         CALL section_vals_get(nonbonded_section, explicit=explicit)
         IF (explicit .AND. qs_control%xtb_control%do_nonbonded) THEN
            CALL section_vals_get(genpot_section, explicit=explicit, n_repetition=ngp)
            IF (explicit) THEN
               CALL pair_potential_reallocate(qs_control%xtb_control%nonbonded, 1, ngp, gp=.true.)
               CALL read_gp_section(qs_control%xtb_control%nonbonded, genpot_section, 0)
            END IF
         END IF !nonbonded
         CALL section_vals_val_get(xtb_section, "EPS_PAIRPOTENTIAL", &
                                   r_val=qs_control%xtb_control%eps_pair)
         ! SR Coulomb
         CALL section_vals_val_get(xtb_parameter, "COULOMB_SR_CUT", r_vals=scal)
         qs_control%xtb_control%coulomb_sr_cut = scal(1)
         CALL section_vals_val_get(xtb_parameter, "COULOMB_SR_EPS", r_vals=scal)
         qs_control%xtb_control%coulomb_sr_eps = scal(1)
         ! XB_radius
         CALL section_vals_val_get(xtb_parameter, "XB_RADIUS", r_val=qs_control%xtb_control%xb_radius)
         ! Kab
         CALL section_vals_val_get(xtb_parameter, "KAB_PARAM", n_rep_val=n_rep)
         ! Coulomb
         SELECT CASE (qs_control%xtb_control%gfn_type)
         CASE (0)
            qs_control%xtb_control%coulomb_interaction = .false.
            qs_control%xtb_control%coulomb_lr = .false.
            qs_control%xtb_control%tb3_interaction = .false.
            qs_control%xtb_control%check_atomic_charges = .false.
            CALL section_vals_val_get(xtb_section, "VARIATIONAL_DIPOLE", &
                                      l_val=qs_control%xtb_control%var_dipole)
         CASE (1)
            ! For debugging purposes
            CALL section_vals_val_get(xtb_section, "COULOMB_INTERACTION", &
                                      l_val=qs_control%xtb_control%coulomb_interaction)
            CALL section_vals_val_get(xtb_section, "COULOMB_LR", &
                                      l_val=qs_control%xtb_control%coulomb_lr)
            CALL section_vals_val_get(xtb_section, "TB3_INTERACTION", &
                                      l_val=qs_control%xtb_control%tb3_interaction)
            ! Check for bad atomic charges
            CALL section_vals_val_get(xtb_section, "CHECK_ATOMIC_CHARGES", &
                                      l_val=qs_control%xtb_control%check_atomic_charges)
            qs_control%xtb_control%var_dipole = .false.
         CASE (2)
            cpabort("gfn2-xtb tbd")
         CASE DEFAULT
            cpabort("GFN type")
         END SELECT
         qs_control%xtb_control%kab_nval = n_rep
         IF (n_rep > 0) THEN
            ALLOCATE (qs_control%xtb_control%kab_param(3, n_rep))
            ALLOCATE (qs_control%xtb_control%kab_types(2, n_rep))
            ALLOCATE (qs_control%xtb_control%kab_vals(n_rep))
            DO j = 1, n_rep
               CALL section_vals_val_get(xtb_parameter, "KAB_PARAM", i_rep_val=j, c_vals=clist)
               qs_control%xtb_control%kab_param(1, j) = clist(1)
               CALL get_ptable_info(clist(1), &
                                    ielement=qs_control%xtb_control%kab_types(1, j))
               qs_control%xtb_control%kab_param(2, j) = clist(2)
               CALL get_ptable_info(clist(2), &
                                    ielement=qs_control%xtb_control%kab_types(2, j))
               qs_control%xtb_control%kab_param(3, j) = clist(3)
               READ (clist(3), '(F10.0)') qs_control%xtb_control%kab_vals(j)
            END DO
         END IF
 
         IF (qs_control%xtb_control%gfn_type == 0) THEN
            CALL section_vals_val_get(xtb_parameter, "SRB_PARAMETER", r_vals=scal)
            qs_control%xtb_control%ksrb = scal(1)
            qs_control%xtb_control%esrb = scal(2)
            qs_control%xtb_control%gscal = scal(3)
            qs_control%xtb_control%c1srb = scal(4)
            qs_control%xtb_control%c2srb = scal(5)
            qs_control%xtb_control%shift = scal(6)
         END IF
 
         CALL section_vals_val_get(xtb_section, "EN_SHIFT_TYPE", c_val=cval)
         CALL uppercase(cval)
         SELECT CASE (trim(cval))
         CASE ("SELECT")
            qs_control%xtb_control%enshift_type = 0
         CASE ("MOLECULE")
            qs_control%xtb_control%enshift_type = 1
         CASE ("CRYSTAL")
            qs_control%xtb_control%enshift_type = 2
         CASE DEFAULT
            cpabort("Unknown value for EN_SHIFT_TYPE")
         END SELECT
 
         ! EEQ solver params
         CALL read_eeq_param(eeq_section, qs_control%xtb_control%eeq_sparam)
      END IF
 
      ! Optimize LRI basis set
      CALL section_vals_get(lri_optbas_section, explicit=qs_control%lri_optbas)
 
      ! Use tblite if selected through XTB/GFN_TYPE TBLITE.
      IF (qs_control%xtb_control%do_tblite) THEN
         CALL section_vals_val_get(xtb_tblite, "METHOD", &
                                   i_val=qs_control%xtb_control%tblite_method)
         CALL section_vals_val_get(xtb_tblite, "PARAM", &
                                   c_val=qs_control%xtb_control%tblite_param_file)
         CALL section_vals_val_get(xtb_tblite, "ACCURACY", &
                                   r_val=qs_control%xtb_control%tblite_accuracy)
         IF (qs_control%xtb_control%tblite_accuracy <= 0.0_dp) &
            cpabort("XTB/TBLITE/ACCURACY must be positive")
         IF (qs_control%xtb_control%tblite_mixer_damping <= 0.0_dp) &
            cpabort("XTB/TBLITE_MIXER/DAMPING must be positive")
         CALL section_vals_val_get(xtb_tblite, "REFERENCE_CLI", l_val=tblite_reference_cli)
         CALL section_vals_get(xtb_tblite_ref_cli, explicit=tblite_reference_cli_section)
         IF (tblite_reference_cli .AND. (.NOT. tblite_reference_cli_section)) &
            cpabort("XTB/TBLITE/REFERENCE_CLI keyword requires an XTB/TBLITE/REFERENCE_CLI section")
         IF (tblite_reference_cli .OR. tblite_reference_cli_section) THEN
            CALL read_xtb_reference_cli_section(xtb_tblite_ref_cli, qs_control%xtb_control%reference_cli, cell)
            qs_control%xtb_control%reference_cli%enabled = .true.
         END IF
         CALL cite_reference(katbashev2025)
         ! tblite handles periodic long-range terms internally from the CP2K cell periodicity.
         ! Keep xtb_control%do_ewald as read above from XTB/DO_EWALD or SUBSYS/CELL/PERIODIC,
         ! matching the DFTB and CP2K-internal xTB setup.
      END IF
 
      CALL timestop(handle)
   END SUBROUTINE read_qs_section
 
! **************************************************************************************************
!> \brief Read a TBLITE_MIXER section.
!> \param mixer_section input section
!> \param iterations tblite SCC iteration limit
!> \param memory Broyden history length
!> \param solver native tblite electronic solver id
!> \param damping mixer damping parameter
!> \param omega0 Broyden regularization weight
!> \param min_weight minimum dynamic Broyden weight
!> \param max_weight maximum dynamic Broyden weight
!> \param weight_factor residual-to-weight scaling factor
!> \param section_name diagnostic section name
! **************************************************************************************************
   SUBROUTINE read_tblite_mixer_section(mixer_section, iterations, memory, solver, damping, omega0, min_weight, &
                                        max_weight, weight_factor, section_name)
      TYPE(section_vals_type), POINTER                   :: mixer_section
      INTEGER, INTENT(INOUT)                             :: iterations, memory, solver
      REAL(kind=dp), INTENT(INOUT)                       :: damping, omega0, min_weight, max_weight, &
                                                            weight_factor
      CHARACTER(LEN=*), INTENT(IN)                       :: section_name
 
      LOGICAL                                            :: explicit, memory_explicit
 
      CALL section_vals_get(mixer_section, explicit=explicit)
      IF (.NOT. explicit) RETURN
 
      CALL section_vals_val_get(mixer_section, "ITERATIONS", i_val=iterations)
      CALL section_vals_val_get(mixer_section, "MEMORY", explicit=memory_explicit)
      IF (memory_explicit) CALL section_vals_val_get(mixer_section, "MEMORY", i_val=memory)
      IF (.NOT. memory_explicit .OR. memory == tblite_mixer_memory_inherit) THEN
         memory = iterations
      END IF
      CALL section_vals_val_get(mixer_section, "SOLVER", i_val=solver)
      CALL section_vals_val_get(mixer_section, "DAMPING", r_val=damping)
      CALL section_vals_val_get(mixer_section, "OMEGA0", r_val=omega0)
      CALL section_vals_val_get(mixer_section, "MIN_WEIGHT", r_val=min_weight)
      CALL section_vals_val_get(mixer_section, "MAX_WEIGHT", r_val=max_weight)
      CALL section_vals_val_get(mixer_section, "WEIGHT_FACTOR", r_val=weight_factor)
 
      IF (iterations < 1) cpabort(trim(section_name)//"/ITERATIONS must be positive")
      IF (memory < 1) cpabort(trim(section_name)//"/MEMORY must be positive")
      SELECT CASE (solver)
      CASE (tblite_solver_gvd, tblite_solver_gvr)
      CASE DEFAULT
         cpabort(trim(section_name)//"/SOLVER must be GVD or GVR")
      END SELECT
      IF (damping <= 0.0_dp) cpabort(trim(section_name)//"/DAMPING must be positive")
      IF (omega0 <= 0.0_dp) cpabort(trim(section_name)//"/OMEGA0 must be positive")
      IF (min_weight <= 0.0_dp) cpabort(trim(section_name)//"/MIN_WEIGHT must be positive")
      IF (max_weight <= 0.0_dp) cpabort(trim(section_name)//"/MAX_WEIGHT must be positive")
      IF (max_weight < min_weight) &
         cpabort(trim(section_name)//"/MAX_WEIGHT must not be smaller than MIN_WEIGHT")
      IF (weight_factor <= 0.0_dp) cpabort(trim(section_name)//"/WEIGHT_FACTOR must be positive")
 
   END SUBROUTINE read_tblite_mixer_section
 
! **************************************************************************************************
!> \brief Read native tblite CLI reference options.
!> \param ref_cli_section input section
!> \param ref_cli reference CLI control data
!> \param cell optional cell used to transform Cartesian input vectors
! **************************************************************************************************
   SUBROUTINE read_xtb_reference_cli_section(ref_cli_section, ref_cli, cell)
      TYPE(section_vals_type), POINTER                   :: ref_cli_section
      TYPE(xtb_reference_cli_type), INTENT(INOUT)        :: ref_cli
      TYPE(cell_type), OPTIONAL, POINTER                 :: cell
 
      REAL(kind=dp), DIMENSION(:), POINTER               :: efield
      TYPE(section_vals_type), POINTER                   :: fit_section, guess_section, &
                                                            param_section, solvation_section, &
                                                            tagdiff_section
 
      CALL section_vals_val_get(ref_cli_section, "_SECTION_PARAMETERS_", l_val=ref_cli%enabled)
      CALL section_vals_val_get(ref_cli_section, "PROGRAM_NAME", c_val=ref_cli%program_name)
      CALL section_vals_val_get(ref_cli_section, "GUESS", i_val=ref_cli%guess)
      CALL section_vals_val_get(ref_cli_section, "WORK_DIRECTORY", c_val=ref_cli%work_directory)
      CALL section_vals_val_get(ref_cli_section, "PREFIX", c_val=ref_cli%prefix)
      CALL section_vals_val_get(ref_cli_section, "INPUT_FORMAT", c_val=ref_cli%input_format)
      CALL section_vals_val_get(ref_cli_section, "RESTART", c_val=ref_cli%restart_file)
      CALL section_vals_val_get(ref_cli_section, "GRAD", c_val=ref_cli%grad_file)
      CALL section_vals_val_get(ref_cli_section, "JSON", c_val=ref_cli%json_file)
      CALL section_vals_val_get(ref_cli_section, "POST_PROCESSING", c_val=ref_cli%post_processing)
      CALL section_vals_val_get(ref_cli_section, "POST_PROCESSING_OUTPUT", &
                                c_val=ref_cli%post_processing_output_file)
      CALL section_vals_val_get(ref_cli_section, "EFIELD", explicit=ref_cli%efield_active)
      IF (ref_cli%efield_active) THEN
         NULLIFY (efield)
         CALL section_vals_val_get(ref_cli_section, "EFIELD", r_vals=efield)
         ref_cli%efield = efield(1:3)
         IF (PRESENT(cell)) THEN
            IF (ASSOCIATED(cell)) CALL cell_transform_input_cartesian(cell, ref_cli%efield)
         END IF
      END IF
      solvation_section => section_vals_get_subs_vals(ref_cli_section, "IMPLICIT_SOLVATION")
      CALL section_vals_get(solvation_section, explicit=ref_cli%solvation_active)
      IF (ref_cli%solvation_active) THEN
         CALL section_vals_val_get(solvation_section, "MODEL", i_val=ref_cli%solvation_model)
         CALL section_vals_val_get(solvation_section, "SOLVENT", c_val=ref_cli%solvation_solvent)
         CALL section_vals_val_get(solvation_section, "BORN_KERNEL", i_val=ref_cli%solvation_born_kernel)
         CALL section_vals_val_get(solvation_section, "SOLUTION_STATE", i_val=ref_cli%solvation_state)
         IF (len_trim(ref_cli%solvation_solvent) == 0) &
            cpabort("REFERENCE_CLI implicit solvation needs SOLVENT")
         IF (ref_cli%solvation_model == tblite_cli_solvation_cpcm .AND. &
             ref_cli%solvation_born_kernel /= tblite_cli_born_kernel_auto) &
            cpabort("BORN_KERNEL is invalid with MODEL CPCM")
         IF (ref_cli%solvation_state /= tblite_cli_solution_state_gsolv) THEN
            SELECT CASE (ref_cli%solvation_model)
            CASE (tblite_cli_solvation_alpb, tblite_cli_solvation_gbsa)
               ! Native tblite supports solution-state shifts for parametrized named-solvent ALPB/GBSA.
            CASE (tblite_cli_solvation_gbe, tblite_cli_solvation_gb, tblite_cli_solvation_cpcm)
               cpabort("SOLUTION_STATE is valid only for ALPB/GBSA")
            END SELECT
         END IF
      END IF
      CALL section_vals_val_get(ref_cli_section, "ELECTRONIC_TEMPERATURE_GUESS", &
                                r_val=ref_cli%electronic_temperature_guess)
      IF (ref_cli%electronic_temperature_guess < 0.0_dp) &
         cpabort("XTB/TBLITE/REFERENCE_CLI/ELECTRONIC_TEMPERATURE_GUESS must not be negative")
      IF (ref_cli%electronic_temperature_guess > 0.0_dp .AND. ref_cli%guess /= tblite_guess_ceh) &
         cpabort("XTB/TBLITE/REFERENCE_CLI/ELECTRONIC_TEMPERATURE_GUESS requires GUESS CEH")
      guess_section => section_vals_get_subs_vals(ref_cli_section, "GUESS_CLI")
      CALL section_vals_get(guess_section, explicit=ref_cli%guess_cli%enabled)
      IF (ref_cli%guess_cli%enabled) THEN
         CALL section_vals_val_get(guess_section, "METHOD", i_val=ref_cli%guess_cli%method)
         CALL section_vals_val_get(guess_section, "ELECTRONIC_TEMPERATURE_GUESS", &
                                   r_val=ref_cli%guess_cli%electronic_temperature_guess)
         IF (ref_cli%guess_cli%electronic_temperature_guess < 0.0_dp) &
            cpabort("REFERENCE_CLI/GUESS_CLI/ELECTRONIC_TEMPERATURE_GUESS must not be negative")
         CALL section_vals_val_get(guess_section, "SOLVER", i_val=ref_cli%guess_cli%solver)
         CALL section_vals_val_get(guess_section, "EFIELD", explicit=ref_cli%guess_cli%efield_active)
         IF (ref_cli%guess_cli%efield_active) THEN
            NULLIFY (efield)
            CALL section_vals_val_get(guess_section, "EFIELD", r_vals=efield)
            ref_cli%guess_cli%efield = efield(1:3)
            IF (PRESENT(cell)) THEN
               IF (ASSOCIATED(cell)) CALL cell_transform_input_cartesian(cell, ref_cli%guess_cli%efield)
            END IF
         END IF
         CALL section_vals_val_get(guess_section, "GRAD", l_val=ref_cli%guess_cli%grad)
         CALL section_vals_val_get(guess_section, "JSON", c_val=ref_cli%guess_cli%json_file)
         CALL section_vals_val_get(guess_section, "INPUT_FORMAT", c_val=ref_cli%guess_cli%input_format)
         CALL section_vals_val_get(guess_section, "INPUT_FILE", c_val=ref_cli%guess_cli%input_file)
      END IF
      param_section => section_vals_get_subs_vals(ref_cli_section, "PARAM_CLI")
      CALL section_vals_get(param_section, explicit=ref_cli%param_cli%enabled)
      IF (ref_cli%param_cli%enabled) THEN
         CALL section_vals_val_get(param_section, "METHOD", explicit=ref_cli%param_cli%method_explicit, &
                                   i_val=ref_cli%param_cli%method)
         CALL section_vals_val_get(param_section, "OUTPUT", c_val=ref_cli%param_cli%output_file)
         CALL section_vals_val_get(param_section, "INPUT_FILE", c_val=ref_cli%param_cli%input_file)
      END IF
      fit_section => section_vals_get_subs_vals(ref_cli_section, "FIT_CLI")
      CALL section_vals_get(fit_section, explicit=ref_cli%fit_cli%enabled)
      IF (ref_cli%fit_cli%enabled) THEN
         CALL section_vals_val_get(fit_section, "PARAM_FILE", c_val=ref_cli%fit_cli%param_file)
         CALL section_vals_val_get(fit_section, "INPUT_FILE", c_val=ref_cli%fit_cli%input_file)
         CALL section_vals_val_get(fit_section, "DRY_RUN", l_val=ref_cli%fit_cli%dry_run)
         CALL section_vals_val_get(fit_section, "COPY", c_val=ref_cli%fit_cli%copy_file)
         IF (len_trim(ref_cli%fit_cli%param_file) == 0) &
            cpabort("XTB/TBLITE/REFERENCE_CLI/FIT_CLI needs PARAM_FILE")
         IF (len_trim(ref_cli%fit_cli%input_file) == 0) &
            cpabort("XTB/TBLITE/REFERENCE_CLI/FIT_CLI needs INPUT_FILE")
      END IF
      tagdiff_section => section_vals_get_subs_vals(ref_cli_section, "TAGDIFF_CLI")
      CALL section_vals_get(tagdiff_section, explicit=ref_cli%tagdiff_cli%enabled)
      IF (ref_cli%tagdiff_cli%enabled) THEN
         CALL section_vals_val_get(tagdiff_section, "ACTUAL", c_val=ref_cli%tagdiff_cli%actual_file)
         CALL section_vals_val_get(tagdiff_section, "REFERENCE", c_val=ref_cli%tagdiff_cli%reference_file)
         CALL section_vals_val_get(tagdiff_section, "FIT", l_val=ref_cli%tagdiff_cli%fit)
         IF (len_trim(ref_cli%tagdiff_cli%actual_file) == 0) &
            cpabort("XTB/TBLITE/REFERENCE_CLI/TAGDIFF_CLI needs ACTUAL")
         IF (len_trim(ref_cli%tagdiff_cli%reference_file) == 0) &
            cpabort("XTB/TBLITE/REFERENCE_CLI/TAGDIFF_CLI needs REFERENCE")
      END IF
      CALL section_vals_val_get(ref_cli_section, "KEEP_FILES", l_val=ref_cli%keep_files)
      CALL section_vals_val_get(ref_cli_section, "ERROR_LIMIT", r_val=ref_cli%error_limit)
      CALL section_vals_val_get(ref_cli_section, "STOP_ON_ERROR", l_val=ref_cli%stop_on_error)
      CALL section_vals_val_get(ref_cli_section, "CHECK_ENERGY", l_val=ref_cli%check_energy)
      CALL section_vals_val_get(ref_cli_section, "CHECK_FORCES", l_val=ref_cli%check_forces)
      CALL section_vals_val_get(ref_cli_section, "CHECK_VIRIAL", l_val=ref_cli%check_virial)
 
   END SUBROUTINE read_xtb_reference_cli_section
 
! **************************************************************************************************
!> \brief Read TDDFPT-related input parameters.
!> \param t_control  TDDFPT control parameters
!> \param t_section  TDDFPT input section
!> \param qs_control Quickstep control parameters
! **************************************************************************************************
   SUBROUTINE read_tddfpt2_control(t_control, t_section, qs_control)
      TYPE(tddfpt2_control_type), POINTER                :: t_control
      TYPE(section_vals_type), POINTER                   :: t_section
      TYPE(qs_control_type), POINTER                     :: qs_control
 
      CHARACTER(LEN=*), PARAMETER :: routinen = 'read_tddfpt2_control'
 
      CHARACTER(LEN=default_string_length), &
         DIMENSION(:), POINTER                           :: tmpstringlist
      INTEGER                                            :: handle, irep, isize, nrep
      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: inds
      LOGICAL                                            :: do_ewald, do_exchange, expl, explicit, &
                                                            multigrid_set
      REAL(kind=dp)                                      :: filter, fval, hfx
      TYPE(section_vals_type), POINTER                   :: dipole_section, mgrid_section, &
                                                            soc_section, stda_section, xc_func, &
                                                            xc_section
 
      CALL timeset(routinen, handle)
 
      CALL section_vals_val_get(t_section, "_SECTION_PARAMETERS_", l_val=t_control%enabled)
 
      CALL section_vals_val_get(t_section, "NSTATES", i_val=t_control%nstates)
      CALL section_vals_val_get(t_section, "MAX_ITER", i_val=t_control%niters)
      CALL section_vals_val_get(t_section, "MAX_KV", i_val=t_control%nkvs)
      CALL section_vals_val_get(t_section, "NLUMO", i_val=t_control%nlumo)
      CALL section_vals_val_get(t_section, "NPROC_STATE", i_val=t_control%nprocs)
      CALL section_vals_val_get(t_section, "KERNEL", i_val=t_control%kernel)
      CALL section_vals_val_get(t_section, "SPINFLIP", i_val=t_control%spinflip)
      CALL section_vals_val_get(t_section, "OE_CORR", i_val=t_control%oe_corr)
      CALL section_vals_val_get(t_section, "EV_SHIFT", r_val=t_control%ev_shift)
      CALL section_vals_val_get(t_section, "EOS_SHIFT", r_val=t_control%eos_shift)
 
      CALL section_vals_val_get(t_section, "CONVERGENCE", r_val=t_control%conv)
      CALL section_vals_val_get(t_section, "MIN_AMPLITUDE", r_val=t_control%min_excitation_amplitude)
      CALL section_vals_val_get(t_section, "ORTHOGONAL_EPS", r_val=t_control%orthogonal_eps)
 
      CALL section_vals_val_get(t_section, "RESTART", l_val=t_control%is_restart)
      CALL section_vals_val_get(t_section, "RKS_TRIPLETS", l_val=t_control%rks_triplets)
      CALL section_vals_val_get(t_section, "DO_LRIGPW", l_val=t_control%do_lrigpw)
      CALL section_vals_val_get(t_section, "DO_SMEARING", l_val=t_control%do_smearing)
      CALL section_vals_val_get(t_section, "DO_BSE", l_val=t_control%do_bse)
      CALL section_vals_val_get(t_section, "DO_BSE_W_ONLY", l_val=t_control%do_bse_w_only)
      CALL section_vals_val_get(t_section, "DO_BSE_GW_ONLY", l_val=t_control%do_bse_gw_only)
      CALL section_vals_val_get(t_section, "ADMM_KERNEL_CORRECTION_SYMMETRIC", l_val=t_control%admm_symm)
      CALL section_vals_val_get(t_section, "ADMM_KERNEL_XC_CORRECTION", l_val=t_control%admm_xc_correction)
      CALL section_vals_val_get(t_section, "EXCITON_DESCRIPTORS", l_val=t_control%do_exciton_descriptors)
      CALL section_vals_val_get(t_section, "DIRECTIONAL_EXCITON_DESCRIPTORS", l_val=t_control%do_directional_exciton_descriptors)
 
      ! read automatically generated auxiliary basis for LRI
      CALL section_vals_val_get(t_section, "AUTO_BASIS", n_rep_val=nrep)
      DO irep = 1, nrep
         CALL section_vals_val_get(t_section, "AUTO_BASIS", i_rep_val=irep, c_vals=tmpstringlist)
         IF (SIZE(tmpstringlist) == 2) THEN
            CALL uppercase(tmpstringlist(2))
            SELECT CASE (tmpstringlist(2))
            CASE ("X")
               isize = -1
            CASE ("SMALL")
               isize = 0
            CASE ("MEDIUM")
               isize = 1
            CASE ("LARGE")
               isize = 2
            CASE ("HUGE")
               isize = 3
            CASE DEFAULT
               cpabort("Unknown basis size in AUTO_BASIS keyword:"//trim(tmpstringlist(1)))
            END SELECT
            !
            SELECT CASE (tmpstringlist(1))
            CASE ("X")
            CASE ("P_LRI_AUX")
               t_control%auto_basis_p_lri_aux = isize
            CASE DEFAULT
               cpabort("Unknown basis type in AUTO_BASIS keyword:"//trim(tmpstringlist(1)))
            END SELECT
         ELSE
            CALL cp_abort(__location__, &
                          "AUTO_BASIS keyword in &PROPERTIES &TDDFT section has a wrong number of arguments.")
         END IF
      END DO
 
      IF (t_control%conv < 0) &
         t_control%conv = abs(t_control%conv)
 
      ! DIPOLE_MOMENTS subsection
      dipole_section => section_vals_get_subs_vals(t_section, "DIPOLE_MOMENTS")
      CALL section_vals_val_get(dipole_section, "DIPOLE_FORM", explicit=explicit)
      IF (explicit) THEN
         CALL section_vals_val_get(dipole_section, "DIPOLE_FORM", i_val=t_control%dipole_form)
      ELSE
         t_control%dipole_form = 0
      END IF
      CALL section_vals_val_get(dipole_section, "REFERENCE", i_val=t_control%dipole_reference)
      CALL section_vals_val_get(dipole_section, "REFERENCE_POINT", explicit=explicit)
      IF (explicit) THEN
         CALL section_vals_val_get(dipole_section, "REFERENCE_POINT", r_vals=t_control%dipole_ref_point)
      ELSE
         NULLIFY (t_control%dipole_ref_point)
         IF (t_control%dipole_form == tddfpt_dipole_length .AND. t_control%dipole_reference == use_mom_ref_user) THEN
            cpabort("User-defined reference point should be given explicitly")
         END IF
      END IF
 
      !SOC subsection
      soc_section => section_vals_get_subs_vals(t_section, "SOC")
      CALL section_vals_get(soc_section, explicit=explicit)
      IF (explicit) THEN
         t_control%do_soc = .true.
      END IF
 
      ! MGRID subsection
      mgrid_section => section_vals_get_subs_vals(t_section, "MGRID")
      CALL section_vals_get(mgrid_section, explicit=t_control%mgrid_is_explicit)
 
      IF (t_control%mgrid_is_explicit) THEN
         CALL section_vals_val_get(mgrid_section, "NGRIDS", i_val=t_control%mgrid_ngrids, explicit=explicit)
         IF (.NOT. explicit) t_control%mgrid_ngrids = SIZE(qs_control%e_cutoff)
 
         CALL section_vals_val_get(mgrid_section, "CUTOFF", r_val=t_control%mgrid_cutoff, explicit=explicit)
         IF (.NOT. explicit) t_control%mgrid_cutoff = qs_control%cutoff
 
         CALL section_vals_val_get(mgrid_section, "PROGRESSION_FACTOR", &
                                   r_val=t_control%mgrid_progression_factor, explicit=explicit)
         IF (explicit) THEN
            IF (t_control%mgrid_progression_factor <= 1.0_dp) &
               CALL cp_abort(__location__, &
                             "Progression factor should be greater then 1.0 to ensure multi-grid ordering")
         ELSE
            t_control%mgrid_progression_factor = qs_control%progression_factor
         END IF
 
         CALL section_vals_val_get(mgrid_section, "COMMENSURATE", l_val=t_control%mgrid_commensurate_mgrids, explicit=explicit)
         IF (.NOT. explicit) t_control%mgrid_commensurate_mgrids = qs_control%commensurate_mgrids
         IF (t_control%mgrid_commensurate_mgrids) THEN
            IF (explicit) THEN
               t_control%mgrid_progression_factor = 4.0_dp
            ELSE
               t_control%mgrid_progression_factor = qs_control%progression_factor
            END IF
         END IF
 
         CALL section_vals_val_get(mgrid_section, "REL_CUTOFF", r_val=t_control%mgrid_relative_cutoff, explicit=explicit)
         IF (.NOT. explicit) t_control%mgrid_relative_cutoff = qs_control%relative_cutoff
 
         CALL section_vals_val_get(mgrid_section, "MULTIGRID_SET", l_val=multigrid_set, explicit=explicit)
         IF (.NOT. explicit) multigrid_set = .false.
         IF (multigrid_set) THEN
            CALL section_vals_val_get(mgrid_section, "MULTIGRID_CUTOFF", r_vals=t_control%mgrid_e_cutoff)
         ELSE
            NULLIFY (t_control%mgrid_e_cutoff)
         END IF
 
         CALL section_vals_val_get(mgrid_section, "REALSPACE", l_val=t_control%mgrid_realspace_mgrids, explicit=explicit)
         IF (.NOT. explicit) t_control%mgrid_realspace_mgrids = qs_control%realspace_mgrids
 
         CALL section_vals_val_get(mgrid_section, "SKIP_LOAD_BALANCE_DISTRIBUTED", &
                                   l_val=t_control%mgrid_skip_load_balance, explicit=explicit)
         IF (.NOT. explicit) t_control%mgrid_skip_load_balance = qs_control%skip_load_balance_distributed
 
         IF (ASSOCIATED(t_control%mgrid_e_cutoff)) THEN
            IF (SIZE(t_control%mgrid_e_cutoff) /= t_control%mgrid_ngrids) &
               cpabort("Inconsistent values for number of multi-grids")
 
            ! sort multi-grids in descending order according to their cutoff values
            t_control%mgrid_e_cutoff = -t_control%mgrid_e_cutoff
            ALLOCATE (inds(t_control%mgrid_ngrids))
            CALL sort(t_control%mgrid_e_cutoff, t_control%mgrid_ngrids, inds)
            DEALLOCATE (inds)
            t_control%mgrid_e_cutoff = -t_control%mgrid_e_cutoff
         END IF
      END IF
 
      ! expand XC subsection (if given explicitly)
      xc_section => section_vals_get_subs_vals(t_section, "XC")
      xc_func => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
      CALL section_vals_get(xc_func, explicit=explicit)
      IF (explicit) &
         CALL xc_functionals_expand(xc_func, xc_section)
 
      ! sTDA subsection
      stda_section => section_vals_get_subs_vals(t_section, "STDA")
      IF (t_control%kernel == tddfpt_kernel_stda) THEN
         t_control%stda_control%hfx_fraction = 0.0_dp
         t_control%stda_control%do_exchange = .true.
         t_control%stda_control%eps_td_filter = 1.e-10_dp
         t_control%stda_control%mn_alpha = -99.0_dp
         t_control%stda_control%mn_beta = -99.0_dp
         ! set default for Ewald method (on/off) dependent on periodicity
         SELECT CASE (qs_control%periodicity)
         CASE (0)
            t_control%stda_control%do_ewald = .false.
         CASE (1)
            t_control%stda_control%do_ewald = .true.
         CASE (2)
            t_control%stda_control%do_ewald = .true.
         CASE (3)
            t_control%stda_control%do_ewald = .true.
         CASE DEFAULT
            cpabort("Illegal value for periodiciy")
         END SELECT
         CALL section_vals_get(stda_section, explicit=explicit)
         IF (explicit) THEN
            CALL section_vals_val_get(stda_section, "HFX_FRACTION", r_val=hfx, explicit=expl)
            IF (expl) t_control%stda_control%hfx_fraction = hfx
            CALL section_vals_val_get(stda_section, "EPS_TD_FILTER", r_val=filter, explicit=expl)
            IF (expl) t_control%stda_control%eps_td_filter = filter
            CALL section_vals_val_get(stda_section, "DO_EWALD", l_val=do_ewald, explicit=expl)
            IF (expl) t_control%stda_control%do_ewald = do_ewald
            CALL section_vals_val_get(stda_section, "DO_EXCHANGE", l_val=do_exchange, explicit=expl)
            IF (expl) t_control%stda_control%do_exchange = do_exchange
            CALL section_vals_val_get(stda_section, "MATAGA_NISHIMOTO_CEXP", r_val=fval)
            t_control%stda_control%mn_alpha = fval
            CALL section_vals_val_get(stda_section, "MATAGA_NISHIMOTO_XEXP", r_val=fval)
            t_control%stda_control%mn_beta = fval
         END IF
         CALL section_vals_val_get(stda_section, "COULOMB_SR_CUT", r_val=fval)
         t_control%stda_control%coulomb_sr_cut = fval
         CALL section_vals_val_get(stda_section, "COULOMB_SR_EPS", r_val=fval)
         t_control%stda_control%coulomb_sr_eps = fval
      END IF
 
      CALL timestop(handle)
   END SUBROUTINE read_tddfpt2_control
 
! **************************************************************************************************
!> \brief Write the DFT control parameters to the output unit.
!> \param dft_control ...
!> \param dft_section ...
! **************************************************************************************************
   SUBROUTINE write_dft_control(dft_control, dft_section)
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(section_vals_type), POINTER                   :: dft_section
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'write_dft_control'
 
      CHARACTER(LEN=20)                                  :: tmpstr
      INTEGER                                            :: handle, i, i_rep, n_rep, output_unit
      REAL(kind=dp)                                      :: density_cut, density_smooth_cut_range, &
                                                            gradient_cut, tau_cut
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(enumeration_type), POINTER                    :: enum
      TYPE(keyword_type), POINTER                        :: keyword
      TYPE(section_type), POINTER                        :: section
      TYPE(section_vals_type), POINTER                   :: xc_section
 
      IF (dft_control%qs_control%semi_empirical) RETURN
      IF (dft_control%qs_control%dftb) RETURN
      IF (dft_control%qs_control%xtb) THEN
         CALL write_xtb_control(dft_control%qs_control%xtb_control, dft_section)
         RETURN
      END IF
      CALL timeset(routinen, handle)
 
      NULLIFY (logger)
      logger => cp_get_default_logger()
 
      output_unit = cp_print_key_unit_nr(logger, dft_section, &
                                         "PRINT%DFT_CONTROL_PARAMETERS", extension=".Log")
 
      IF (output_unit > 0) THEN
 
         xc_section => section_vals_get_subs_vals(dft_section, "XC")
 
         IF (dft_control%uks) THEN
            WRITE (unit=output_unit, fmt="(/,T2,A,T78,A)") &
               "DFT| Spin unrestricted (spin-polarized) Kohn-Sham calculation", "UKS"
         ELSE IF (dft_control%roks) THEN
            WRITE (unit=output_unit, fmt="(/,T2,A,T77,A)") &
               "DFT| Spin restricted open Kohn-Sham calculation", "ROKS"
         ELSE
            WRITE (unit=output_unit, fmt="(/,T2,A,T78,A)") &
               "DFT| Spin restricted Kohn-Sham (RKS) calculation", "RKS"
         END IF
 
         WRITE (unit=output_unit, fmt="(T2,A,T76,I5)") &
            "DFT| Multiplicity", dft_control%multiplicity
         WRITE (unit=output_unit, fmt="(T2,A,T76,I5)") &
            "DFT| Number of spin states", dft_control%nspins
 
         WRITE (unit=output_unit, fmt="(T2,A,T76,I5)") &
            "DFT| Charge", dft_control%charge
 
         IF (dft_control%sic_method_id /= sic_none) CALL cite_reference(vandevondele2005b)
         SELECT CASE (dft_control%sic_method_id)
         CASE (sic_none)
            tmpstr = "NO"
         CASE (sic_mauri_spz)
            tmpstr = "SPZ/MAURI SIC"
         CASE (sic_mauri_us)
            tmpstr = "US/MAURI SIC"
         CASE (sic_ad)
            tmpstr = "AD SIC"
         CASE (sic_eo)
            tmpstr = "Explicit Orbital SIC"
         CASE DEFAULT
            ! fix throughout the cp2k for this option
            cpabort("SIC option unknown")
         END SELECT
 
         WRITE (unit=output_unit, fmt="(T2,A,T61,A20)") &
            "DFT| Self-interaction correction (SIC)", adjustr(trim(tmpstr))
 
         IF (dft_control%sic_method_id /= sic_none) THEN
            WRITE (unit=output_unit, fmt="(T2,A,T66,ES15.6)") &
               "DFT| SIC scaling parameter a", dft_control%sic_scaling_a, &
               "DFT| SIC scaling parameter b", dft_control%sic_scaling_b
         END IF
 
         IF (dft_control%sic_method_id == sic_eo) THEN
            IF (dft_control%sic_list_id == sic_list_all) THEN
               WRITE (unit=output_unit, fmt="(T2,A,T66,A)") &
                  "DFT| SIC orbitals", "ALL"
            END IF
            IF (dft_control%sic_list_id == sic_list_unpaired) THEN
               WRITE (unit=output_unit, fmt="(T2,A,T66,A)") &
                  "DFT| SIC orbitals", "UNPAIRED"
            END IF
         END IF
 
         CALL section_vals_val_get(xc_section, "density_cutoff", r_val=density_cut)
         CALL section_vals_val_get(xc_section, "gradient_cutoff", r_val=gradient_cut)
         CALL section_vals_val_get(xc_section, "tau_cutoff", r_val=tau_cut)
         CALL section_vals_val_get(xc_section, "density_smooth_cutoff_range", r_val=density_smooth_cut_range)
 
         WRITE (unit=output_unit, fmt="(T2,A,T66,ES15.6)") &
            "DFT| Cutoffs: density ", density_cut, &
            "DFT|          gradient", gradient_cut, &
            "DFT|          tau     ", tau_cut, &
            "DFT|          cutoff_smoothing_range", density_smooth_cut_range
         CALL section_vals_val_get(xc_section, "XC_GRID%XC_SMOOTH_RHO", &
                                   c_val=tmpstr)
         WRITE (output_unit, '( A, T61, A )') &
            " DFT| XC density smoothing ", adjustr(tmpstr)
         CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
                                   c_val=tmpstr)
         WRITE (output_unit, '( A, T61, A )') &
            " DFT| XC derivatives ", adjustr(tmpstr)
         IF (dft_control%dft_plus_u) THEN
            NULLIFY (enum, keyword, section)
            CALL create_dft_section(section)
            keyword => section_get_keyword(section, "PLUS_U_METHOD")
            CALL keyword_get(keyword, enum=enum)
            WRITE (unit=output_unit, fmt="(/,T2,A,T41,A40)") &
               "DFT+U| Method", adjustr(trim(enum_i2c(enum, dft_control%plus_u_method_id)))
            WRITE (unit=output_unit, fmt="(T2,A)") &
               "DFT+U| Check atomic kind information for details"
            CALL section_release(section)
         END IF
 
         WRITE (unit=output_unit, fmt="(A)") ""
         CALL xc_write(output_unit, xc_section, dft_control%lsd)
 
         IF (dft_control%apply_period_efield) THEN
            WRITE (unit=output_unit, fmt="(A)") ""
            IF (dft_control%period_efield%displacement_field) THEN
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "PERIODIC_EFIELD| Use displacement field formulation"
               WRITE (unit=output_unit, fmt="(T2,A,T66,1X,ES14.6)") &
                  "PERIODIC_EFIELD| Displacement field filter: x", &
                  dft_control%period_efield%d_filter(1), &
                  "PERIODIC_EFIELD|                            y", &
                  dft_control%period_efield%d_filter(2), &
                  "PERIODIC_EFIELD|                            z", &
                  dft_control%period_efield%d_filter(3)
            END IF
            WRITE (unit=output_unit, fmt="(T2,A,T66,1X,ES14.6)") &
               "PERIODIC_EFIELD| Polarisation vector:       x", &
               dft_control%period_efield%polarisation(1), &
               "PERIODIC_EFIELD|                            y", &
               dft_control%period_efield%polarisation(2), &
               "PERIODIC_EFIELD|                            z", &
               dft_control%period_efield%polarisation(3)
 
            WRITE (unit=output_unit, fmt="(T2,A,T66,1X,I14)") &
               "PERIODIC_EFIELD| Start Frame:", &
               dft_control%period_efield%start_frame, &
               "PERIODIC_EFIELD| End Frame:", &
               dft_control%period_efield%end_frame
 
            IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
               WRITE (unit=output_unit, fmt="(T2,A,T66,1X,I14)") &
                  "PERIODIC_EFIELD| Number of Intensities:", &
                  SIZE(dft_control%period_efield%strength_list)
               WRITE (unit=output_unit, fmt="(T2,A,I10,T66,1X,ES14.6)") &
                  "PERIODIC_EFIELD| Intensity List [a.u.] ", &
                  1, dft_control%period_efield%strength_list(1)
               DO i = 2, SIZE(dft_control%period_efield%strength_list)
                  WRITE (unit=output_unit, fmt="(T2,A,I10,T66,1X,ES14.6)") &
                     "PERIODIC_EFIELD|                       ", &
                     i, dft_control%period_efield%strength_list(i)
               END DO
            ELSE
               WRITE (unit=output_unit, fmt="(T2,A,T66,1X,ES14.6)") &
                  "PERIODIC_EFIELD| Intensity [a.u.]:", &
                  dft_control%period_efield%strength
            END IF
 
            IF (sqrt(dot_product(dft_control%period_efield%polarisation, &
                                 dft_control%period_efield%polarisation)) < epsilon(0.0_dp)) THEN
               cpabort("Invalid (too small) polarisation vector specified for PERIODIC_EFIELD")
            END IF
         END IF
 
         IF (dft_control%do_sccs) THEN
            WRITE (unit=output_unit, fmt="(/,T2,A)") &
               "SCCS| Self-consistent continuum solvation model"
            WRITE (unit=output_unit, fmt="(T2,A,T61,ES20.6)") &
               "SCCS| Relative permittivity of the solvent (medium)", &
               dft_control%sccs_control%epsilon_solvent, &
               "SCCS| Absolute permittivity [a.u.]", &
               dft_control%sccs_control%epsilon_solvent/fourpi
            SELECT CASE (dft_control%sccs_control%method_id)
            CASE (sccs_andreussi)
               WRITE (unit=output_unit, fmt="(T2,A,/,(T2,A,T61,ES20.6))") &
                  "SCCS| Dielectric function proposed by Andreussi et al.", &
                  "SCCS|  rho_max", dft_control%sccs_control%rho_max, &
                  "SCCS|  rho_min", dft_control%sccs_control%rho_min
            CASE (sccs_fattebert_gygi)
               WRITE (unit=output_unit, fmt="(T2,A,/,(T2,A,T61,ES20.6))") &
                  "SCCS| Dielectric function proposed by Fattebert and Gygi", &
                  "SCCS|  beta", dft_control%sccs_control%beta, &
                  "SCCS|  rho_zero", dft_control%sccs_control%rho_zero
            CASE DEFAULT
               cpabort("Invalid SCCS model specified. Please, check your input!")
            END SELECT
            SELECT CASE (dft_control%sccs_control%derivative_method)
            CASE (sccs_derivative_fft)
               WRITE (unit=output_unit, fmt="(T2,A,T46,A35)") &
                  "SCCS| Numerical derivative calculation", &
                  adjustr("FFT")
            CASE (sccs_derivative_cd3)
               WRITE (unit=output_unit, fmt="(T2,A,T46,A35)") &
                  "SCCS| Numerical derivative calculation", &
                  adjustr("3-point stencil central differences")
            CASE (sccs_derivative_cd5)
               WRITE (unit=output_unit, fmt="(T2,A,T46,A35)") &
                  "SCCS| Numerical derivative calculation", &
                  adjustr("5-point stencil central differences")
            CASE (sccs_derivative_cd7)
               WRITE (unit=output_unit, fmt="(T2,A,T46,A35)") &
                  "SCCS| Numerical derivative calculation", &
                  adjustr("7-point stencil central differences")
            CASE DEFAULT
               CALL cp_abort(__location__, &
                             "Invalid derivative method specified for SCCS model. "// &
                             "Please, check your input!")
            END SELECT
            WRITE (unit=output_unit, fmt="(T2,A,T61,ES20.6)") &
               "SCCS| Repulsion parameter alpha [mN/m] = [dyn/cm]", &
               cp_unit_from_cp2k(dft_control%sccs_control%alpha_solvent, "mN/m")
            WRITE (unit=output_unit, fmt="(T2,A,T61,ES20.6)") &
               "SCCS| Dispersion parameter beta [GPa]", &
               cp_unit_from_cp2k(dft_control%sccs_control%beta_solvent, "GPa")
            WRITE (unit=output_unit, fmt="(T2,A,T61,ES20.6)") &
               "SCCS| Surface tension gamma [mN/m] = [dyn/cm]", &
               cp_unit_from_cp2k(dft_control%sccs_control%gamma_solvent, "mN/m")
            WRITE (unit=output_unit, fmt="(T2,A,T61,ES20.6)") &
               "SCCS| Mixing parameter applied during the iteration cycle", &
               dft_control%sccs_control%mixing
            WRITE (unit=output_unit, fmt="(T2,A,T61,ES20.6)") &
               "SCCS| Tolerance for the convergence of the SCCS iteration cycle", &
               dft_control%sccs_control%eps_sccs
            WRITE (unit=output_unit, fmt="(T2,A,T61,I20)") &
               "SCCS| Maximum number of iteration steps", &
               dft_control%sccs_control%max_iter
            WRITE (unit=output_unit, fmt="(T2,A,T61,ES20.6)") &
               "SCCS| SCF convergence threshold for starting the SCCS iteration", &
               dft_control%sccs_control%eps_scf
            WRITE (unit=output_unit, fmt="(T2,A,T61,ES20.6)") &
               "SCCS| Numerical increment for the cavity surface calculation", &
               dft_control%sccs_control%delta_rho
         END IF
 
         WRITE (unit=output_unit, fmt="(A)") ""
 
      END IF
 
      IF (dft_control%hairy_probes .EQV. .true.) THEN
         n_rep = SIZE(dft_control%probe)
         IF (output_unit > 0) THEN
            DO i_rep = 1, n_rep
               WRITE (unit=output_unit, fmt="(T2,A,I5)") &
                  "HP | hair probe set", i_rep
               WRITE (unit=output_unit, fmt="(T2,A,T61,*(I5))") &
                  "HP| atom indexes", &
                  (dft_control%probe(i_rep)%atom_ids(i), i=1, dft_control%probe(i_rep)%natoms)
               WRITE (unit=output_unit, fmt="(T2,A,T61,ES20.6)") &
                  "HP| potential", dft_control%probe(i_rep)%mu
               WRITE (unit=output_unit, fmt="(T2,A,T61,F20.2)") &
                  "HP| temperature", dft_control%probe(i_rep)%T
               WRITE (unit=output_unit, fmt="(T2,A,T61,ES20.6)") &
                  "HP| eps_hp", dft_control%probe(i_rep)%eps_hp
            END DO
         END IF
      END IF
 
      CALL cp_print_key_finished_output(output_unit, logger, dft_section, &
                                        "PRINT%DFT_CONTROL_PARAMETERS")
 
      CALL timestop(handle)
 
   END SUBROUTINE write_dft_control
 
! **************************************************************************************************
!> \brief Write the ADMM control parameters to the output unit.
!> \param admm_control ...
!> \param dft_section ...
! **************************************************************************************************
   SUBROUTINE write_admm_control(admm_control, dft_section)
      TYPE(admm_control_type), POINTER                   :: admm_control
      TYPE(section_vals_type), POINTER                   :: dft_section
 
      INTEGER                                            :: iounit
      TYPE(cp_logger_type), POINTER                      :: logger
 
      NULLIFY (logger)
      logger => cp_get_default_logger()
 
      iounit = cp_print_key_unit_nr(logger, dft_section, &
                                    "PRINT%DFT_CONTROL_PARAMETERS", extension=".Log")
 
      IF (iounit > 0) THEN
 
         SELECT CASE (admm_control%admm_type)
         CASE (no_admm_type)
            WRITE (unit=iounit, fmt="(/,T2,A,T77,A)") "ADMM| Specific ADMM type specified", "NONE"
         CASE (admm1_type)
            WRITE (unit=iounit, fmt="(/,T2,A,T76,A)") "ADMM| Specific ADMM type specified", "ADMM1"
         CASE (admm2_type)
            WRITE (unit=iounit, fmt="(/,T2,A,T76,A)") "ADMM| Specific ADMM type specified", "ADMM2"
         CASE (admms_type)
            WRITE (unit=iounit, fmt="(/,T2,A,T76,A)") "ADMM| Specific ADMM type specified", "ADMMS"
         CASE (admmp_type)
            WRITE (unit=iounit, fmt="(/,T2,A,T76,A)") "ADMM| Specific ADMM type specified", "ADMMP"
         CASE (admmq_type)
            WRITE (unit=iounit, fmt="(/,T2,A,T76,A)") "ADMM| Specific ADMM type specified", "ADMMQ"
         CASE DEFAULT
            cpabort("admm_type")
         END SELECT
 
         SELECT CASE (admm_control%purification_method)
         CASE (do_admm_purify_none)
            WRITE (unit=iounit, fmt="(T2,A,T77,A)") "ADMM| Density matrix purification method", "NONE"
         CASE (do_admm_purify_cauchy)
            WRITE (unit=iounit, fmt="(T2,A,T75,A)") "ADMM| Density matrix purification method", "Cauchy"
         CASE (do_admm_purify_cauchy_subspace)
            WRITE (unit=iounit, fmt="(T2,A,T66,A)") "ADMM| Density matrix purification method", "Cauchy subspace"
         CASE (do_admm_purify_mo_diag)
            WRITE (unit=iounit, fmt="(T2,A,T63,A)") "ADMM| Density matrix purification method", "MO diagonalization"
         CASE (do_admm_purify_mo_no_diag)
            WRITE (unit=iounit, fmt="(T2,A,T71,A)") "ADMM| Density matrix purification method", "MO no diag"
         CASE (do_admm_purify_mcweeny)
            WRITE (unit=iounit, fmt="(T2,A,T74,A)") "ADMM| Density matrix purification method", "McWeeny"
         CASE (do_admm_purify_none_dm)
            WRITE (unit=iounit, fmt="(T2,A,T73,A)") "ADMM| Density matrix purification method", "NONE(DM)"
         CASE DEFAULT
            cpabort("admm_purification_method")
         END SELECT
 
         SELECT CASE (admm_control%method)
         CASE (do_admm_basis_projection)
            WRITE (unit=iounit, fmt="(T2,A)") "ADMM| Orbital projection on ADMM basis"
         CASE (do_admm_blocking_purify_full)
            WRITE (unit=iounit, fmt="(T2,A)") "ADMM| Blocked Fock matrix projection with full purification"
         CASE (do_admm_blocked_projection)
            WRITE (unit=iounit, fmt="(T2,A)") "ADMM| Blocked Fock matrix projection"
         CASE (do_admm_charge_constrained_projection)
            WRITE (unit=iounit, fmt="(T2,A)") "ADMM| Orbital projection with charge constrain"
         CASE DEFAULT
            cpabort("admm method")
         END SELECT
 
         SELECT CASE (admm_control%scaling_model)
         CASE (do_admm_exch_scaling_none)
         CASE (do_admm_exch_scaling_merlot)
            WRITE (unit=iounit, fmt="(T2,A)") "ADMM| Use Merlot (2014) scaling model"
         CASE DEFAULT
            cpabort("admm scaling_model")
         END SELECT
 
         WRITE (unit=iounit, fmt="(T2,A,T61,G20.10)") "ADMM| eps_filter", admm_control%eps_filter
 
         SELECT CASE (admm_control%aux_exch_func)
         CASE (do_admm_aux_exch_func_none)
            WRITE (unit=iounit, fmt="(T2,A)") "ADMM| No exchange functional correction term used"
         CASE (do_admm_aux_exch_func_default, do_admm_aux_exch_func_default_libxc)
            WRITE (unit=iounit, fmt="(T2,A,T74,A)") "ADMM| Exchange functional in correction term", "(W)PBEX"
         CASE (do_admm_aux_exch_func_pbex, do_admm_aux_exch_func_pbex_libxc)
            WRITE (unit=iounit, fmt="(T2,A,T77,A)") "ADMM| Exchange functional in correction term", "PBEX"
         CASE (do_admm_aux_exch_func_opt, do_admm_aux_exch_func_opt_libxc)
            WRITE (unit=iounit, fmt="(T2,A,T77,A)") "ADMM| Exchange functional in correction term", "OPTX"
         CASE (do_admm_aux_exch_func_bee, do_admm_aux_exch_func_bee_libxc)
            WRITE (unit=iounit, fmt="(T2,A,T74,A)") "ADMM| Exchange functional in correction term", "Becke88"
         CASE (do_admm_aux_exch_func_sx_libxc)
            WRITE (unit=iounit, fmt="(T2,A,T74,A)") "ADMM| Exchange functional in correction term", "SlaterX"
         CASE DEFAULT
            cpabort("admm aux_exch_func")
         END SELECT
 
         WRITE (unit=iounit, fmt="(A)") ""
 
      END IF
 
      CALL cp_print_key_finished_output(iounit, logger, dft_section, &
                                        "PRINT%DFT_CONTROL_PARAMETERS")
   END SUBROUTINE write_admm_control
 
! **************************************************************************************************
!> \brief Write the xTB control parameters to the output unit.
!> \param xtb_control ...
!> \param dft_section ...
! **************************************************************************************************
   SUBROUTINE write_xtb_control(xtb_control, dft_section)
      TYPE(xtb_control_type), POINTER                    :: xtb_control
      TYPE(section_vals_type), POINTER                   :: dft_section
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'write_xtb_control'
 
      CHARACTER(LEN=16)                                  :: scc_mixer_name, solver_name
      INTEGER                                            :: handle, output_unit
      TYPE(cp_logger_type), POINTER                      :: logger
 
      CALL timeset(routinen, handle)
      NULLIFY (logger)
      logger => cp_get_default_logger()
 
      output_unit = cp_print_key_unit_nr(logger, dft_section, &
                                         "PRINT%DFT_CONTROL_PARAMETERS", extension=".Log")
 
      IF (output_unit > 0) THEN
 
         WRITE (unit=output_unit, fmt="(/,T2,A,T31,A50)") &
            "xTB| Parameter file", adjustr(trim(xtb_control%parameter_file_name))
         WRITE (unit=output_unit, fmt="(T2,A,T71,I10)") &
            "xTB| Basis expansion STO-NG", xtb_control%sto_ng
         WRITE (unit=output_unit, fmt="(T2,A,T71,I10)") &
            "xTB| Basis expansion STO-NG for Hydrogen", xtb_control%h_sto_ng
         WRITE (unit=output_unit, fmt="(T2,A,T71,E10.4)") &
            "xTB| Repulsive pair potential accuracy", xtb_control%eps_pair
         WRITE (unit=output_unit, fmt="(T2,A,T71,F10.6)") &
            "xTB| Repulsive enhancement factor", xtb_control%enscale
         WRITE (unit=output_unit, fmt="(T2,A,T71,L10)") &
            "xTB| Halogen interaction potential", xtb_control%xb_interaction
         WRITE (unit=output_unit, fmt="(T2,A,T71,F10.3)") &
            "xTB| Halogen interaction potential cutoff radius", xtb_control%xb_radius
         WRITE (unit=output_unit, fmt="(T2,A,T71,L10)") &
            "xTB| Nonbonded interactions", xtb_control%do_nonbonded
         SELECT CASE (xtb_control%vdw_type)
         CASE (xtb_vdw_type_none)
            WRITE (unit=output_unit, fmt="(T2,A)") "xTB| No vdW potential selected"
         CASE (xtb_vdw_type_d3)
            WRITE (unit=output_unit, fmt="(T2,A,T72,A)") "xTB| vdW potential type:", "DFTD3(BJ)"
            WRITE (unit=output_unit, fmt="(T2,A,T31,A50)") &
               "xTB| D3 Dispersion: Parameter file", adjustr(trim(xtb_control%dispersion_parameter_file))
         CASE (xtb_vdw_type_d4)
            WRITE (unit=output_unit, fmt="(T2,A,T76,A)") "xTB| vdW potential type:", "DFTD4"
            WRITE (unit=output_unit, fmt="(T2,A,T31,A50)") &
               "xTB| D4 Dispersion: Parameter file", adjustr(trim(xtb_control%dispersion_parameter_file))
         CASE DEFAULT
            cpabort("vdw type")
         END SELECT
         WRITE (unit=output_unit, fmt="(T2,A,T51,3F10.3)") &
            "xTB| Huckel constants ks kp kd", xtb_control%ks, xtb_control%kp, xtb_control%kd
         WRITE (unit=output_unit, fmt="(T2,A,T61,2F10.3)") &
            "xTB| Huckel constants ksp k2sh", xtb_control%ksp, xtb_control%k2sh
         WRITE (unit=output_unit, fmt="(T2,A,T71,F10.3)") &
            "xTB| Mataga-Nishimoto exponent", xtb_control%kg
         WRITE (unit=output_unit, fmt="(T2,A,T71,F10.3)") &
            "xTB| Repulsion potential exponent", xtb_control%kf
         WRITE (unit=output_unit, fmt="(T2,A,T51,3F10.3)") &
            "xTB| Coordination number scaling kcn(s) kcn(p) kcn(d)", &
            xtb_control%kcns, xtb_control%kcnp, xtb_control%kcnd
         WRITE (unit=output_unit, fmt="(T2,A,T71,F10.3)") &
            "xTB| Electronegativity scaling", xtb_control%ken
         WRITE (unit=output_unit, fmt="(T2,A,T61,2F10.3)") &
            "xTB| Halogen potential scaling kxr kx2", xtb_control%kxr, xtb_control%kx2
         SELECT CASE (xtb_control%tblite_scc_mixer)
         CASE (tblite_scc_mixer_auto)
            scc_mixer_name = "AUTO"
         CASE (tblite_scc_mixer_tblite)
            scc_mixer_name = "TBLITE"
         CASE (tblite_scc_mixer_cp2k)
            scc_mixer_name = "CP2K"
         CASE (tblite_scc_mixer_none)
            scc_mixer_name = "NONE"
         CASE DEFAULT
            cpabort("Unknown tblite SCC mixer")
         END SELECT
         SELECT CASE (xtb_control%tblite_mixer_solver)
         CASE (tblite_solver_gvd)
            solver_name = "GVD"
         CASE (tblite_solver_gvr)
            solver_name = "GVR"
         CASE DEFAULT
            cpabort("Unknown tblite SCC mixer solver")
         END SELECT
         WRITE (unit=output_unit, fmt="(T2,A,T72,A)") &
            "xTB| SCC mixer:", trim(scc_mixer_name)
         WRITE (unit=output_unit, fmt="(T2,A,T71,ES10.3)") &
            "xTB| tblite accuracy:", xtb_control%tblite_accuracy
         IF (len_trim(xtb_control%tblite_param_file) > 0) THEN
            WRITE (unit=output_unit, fmt="(T2,A,T33,A)") &
               "xTB| tblite parameter file:", trim(xtb_control%tblite_param_file)
         END IF
         WRITE (unit=output_unit, fmt="(T2,A,T71,F10.3)") &
            "xTB| tblite SCC mixer damping:", xtb_control%tblite_mixer_damping
         WRITE (unit=output_unit, fmt="(T2,A,T71,I10)") &
            "xTB| tblite SCC mixer iterations:", xtb_control%tblite_mixer_iterations
         WRITE (unit=output_unit, fmt="(T2,A,T71,I10)") &
            "xTB| tblite SCC mixer memory:", xtb_control%tblite_mixer_memory
         WRITE (unit=output_unit, fmt="(T2,A,T72,A)") &
            "xTB| tblite SCC mixer solver:", trim(solver_name)
         WRITE (unit=output_unit, fmt="(T2,A,T71,ES10.3)") &
            "xTB| tblite SCC mixer omega0:", xtb_control%tblite_mixer_omega0
         WRITE (unit=output_unit, fmt="(T2,A,T71,ES10.3)") &
            "xTB| tblite SCC mixer min weight:", xtb_control%tblite_mixer_min_weight
         WRITE (unit=output_unit, fmt="(T2,A,T71,ES10.3)") &
            "xTB| tblite SCC mixer max weight:", xtb_control%tblite_mixer_max_weight
         WRITE (unit=output_unit, fmt="(T2,A,T71,ES10.3)") &
            "xTB| tblite SCC mixer weight factor:", xtb_control%tblite_mixer_weight_factor
         WRITE (unit=output_unit, fmt="(/)")
 
      END IF
 
      CALL cp_print_key_finished_output(output_unit, logger, dft_section, &
                                        "PRINT%DFT_CONTROL_PARAMETERS")
 
      CALL timestop(handle)
 
   END SUBROUTINE write_xtb_control
 
! **************************************************************************************************
!> \brief Purpose: Write the QS control parameters to the output unit.
!> \param qs_control ...
!> \param dft_section ...
! **************************************************************************************************
   SUBROUTINE write_qs_control(qs_control, dft_section)
      TYPE(qs_control_type), INTENT(IN)                  :: qs_control
      TYPE(section_vals_type), POINTER                   :: dft_section
 
      CHARACTER(len=*), PARAMETER                        :: routinen = 'write_qs_control'
 
      CHARACTER(len=20)                                  :: method, quadrature
      INTEGER                                            :: handle, i, igrid_level, ngrid_level, &
                                                            output_unit
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(ddapc_restraint_type), POINTER                :: ddapc_restraint_control
      TYPE(enumeration_type), POINTER                    :: enum
      TYPE(keyword_type), POINTER                        :: keyword
      TYPE(section_type), POINTER                        :: qs_section
      TYPE(section_vals_type), POINTER                   :: print_section_vals, qs_section_vals
 
      IF (qs_control%semi_empirical) RETURN
      IF (qs_control%dftb) RETURN
      IF (qs_control%xtb) RETURN
      CALL timeset(routinen, handle)
      NULLIFY (logger, print_section_vals, qs_section, qs_section_vals)
      logger => cp_get_default_logger()
      print_section_vals => section_vals_get_subs_vals(dft_section, "PRINT")
      qs_section_vals => section_vals_get_subs_vals(dft_section, "QS")
      CALL section_vals_get(qs_section_vals, section=qs_section)
 
      NULLIFY (enum, keyword)
      keyword => section_get_keyword(qs_section, "METHOD")
      CALL keyword_get(keyword, enum=enum)
      method = trim(enum_i2c(enum, qs_control%method_id))
 
      NULLIFY (enum, keyword)
      keyword => section_get_keyword(qs_section, "QUADRATURE")
      CALL keyword_get(keyword, enum=enum)
      quadrature = trim(enum_i2c(enum, qs_control%gapw_control%quadrature))
 
      output_unit = cp_print_key_unit_nr(logger, print_section_vals, &
                                         "DFT_CONTROL_PARAMETERS", extension=".Log")
      IF (output_unit > 0) THEN
         ngrid_level = SIZE(qs_control%e_cutoff)
         WRITE (unit=output_unit, fmt="(/,T2,A,T61,A20)") &
            "QS| Method:", adjustr(method)
         IF (qs_control%pw_grid_opt%spherical) THEN
            WRITE (unit=output_unit, fmt="(T2,A,T61,A)") &
               "QS| Density plane wave grid type", " SPHERICAL HALFSPACE"
         ELSE IF (qs_control%pw_grid_opt%fullspace) THEN
            WRITE (unit=output_unit, fmt="(T2,A,T57,A)") &
               "QS| Density plane wave grid type", " NON-SPHERICAL FULLSPACE"
         ELSE
            WRITE (unit=output_unit, fmt="(T2,A,T57,A)") &
               "QS| Density plane wave grid type", " NON-SPHERICAL HALFSPACE"
         END IF
         WRITE (unit=output_unit, fmt="(T2,A,T71,I10)") &
            "QS| Number of grid levels:", SIZE(qs_control%e_cutoff)
         IF (ngrid_level == 1) THEN
            WRITE (unit=output_unit, fmt="(T2,A,T71,F10.1)") &
               "QS| Density cutoff [a.u.]:", qs_control%e_cutoff(1)
         ELSE
            WRITE (unit=output_unit, fmt="(T2,A,T71,F10.1)") &
               "QS| Density cutoff [a.u.]:", qs_control%cutoff
            IF (qs_control%commensurate_mgrids) &
               WRITE (unit=output_unit, fmt="(T2,A)") "QS| Using commensurate multigrids"
            WRITE (unit=output_unit, fmt="(T2,A,T71,F10.1)") &
               "QS| Multi grid cutoff [a.u.]: 1) grid level", qs_control%e_cutoff(1)
            WRITE (unit=output_unit, fmt="(T2,A,I3,A,T71,F10.1)") &
               ("QS|                         ", igrid_level, ") grid level", &
                qs_control%e_cutoff(igrid_level), &
                igrid_level=2, SIZE(qs_control%e_cutoff))
         END IF
         IF (qs_control%pao) THEN
            WRITE (unit=output_unit, fmt="(T2,A)") "QS| PAO active"
         END IF
         WRITE (unit=output_unit, fmt="(T2,A,T71,F10.1)") &
            "QS| Grid level progression factor:", qs_control%progression_factor
         WRITE (unit=output_unit, fmt="(T2,A,T71,F10.1)") &
            "QS| Relative density cutoff [a.u.]:", qs_control%relative_cutoff
         WRITE (unit=output_unit, fmt="(T2,A,T73,ES8.1)") &
            "QS| Interaction thresholds: eps_pgf_orb:", &
            qs_control%eps_pgf_orb, &
            "QS|                         eps_filter_matrix:", &
            qs_control%eps_filter_matrix, &
            "QS|                         eps_core_charge:", &
            qs_control%eps_core_charge, &
            "QS|                         eps_rho_gspace:", &
            qs_control%eps_rho_gspace, &
            "QS|                         eps_rho_rspace:", &
            qs_control%eps_rho_rspace, &
            "QS|                         eps_gvg_rspace:", &
            qs_control%eps_gvg_rspace, &
            "QS|                         eps_ppl:", &
            qs_control%eps_ppl, &
            "QS|                         eps_ppnl:", &
            qs_control%eps_ppnl
         IF (qs_control%gapw) THEN
            IF (qs_control%gapw_control%accurate_xcint) THEN
               WRITE (unit=output_unit, fmt="(T2,A,T69,F12.6)") &
                  "QS| GAPW|      XC integration using accurate scheme: Ref. exponent =", &
                  qs_control%gapw_control%aweights
            END IF
            !
            SELECT CASE (qs_control%gapw_control%basis_1c)
            CASE (gapw_1c_orb)
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "QS| GAPW|      One center basis from orbital basis primitives"
            CASE (gapw_1c_small)
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "QS| GAPW|      One center basis extended with primitives (small:s)"
            CASE (gapw_1c_medium)
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "QS| GAPW|      One center basis extended with primitives (medium:sp)"
            CASE (gapw_1c_large)
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "QS| GAPW|      One center basis extended with primitives (large:spd)"
            CASE (gapw_1c_very_large)
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "QS| GAPW|      One center basis extended with primitives (very large:spdf)"
            CASE DEFAULT
               cpabort("basis_1c incorrect")
            END SELECT
            WRITE (unit=output_unit, fmt="(T2,A,T73,ES8.1)") &
               "QS| GAPW|                   eps_fit:", &
               qs_control%gapw_control%eps_fit, &
               "QS| GAPW|                   eps_iso:", &
               qs_control%gapw_control%eps_iso, &
               "QS| GAPW|                   eps_svd:", &
               qs_control%gapw_control%eps_svd, &
               "QS| GAPW|                   eps_cpc:", &
               qs_control%gapw_control%eps_cpc
            WRITE (unit=output_unit, fmt="(T2,A,T61,A20)") &
               "QS| GAPW|   atom-r-grid: quadrature:", &
               adjustr(quadrature)
            WRITE (unit=output_unit, fmt="(T2,A,T71,I10)") &
               "QS| GAPW|      atom-s-grid:  max l :", &
               qs_control%gapw_control%lmax_sphere, &
               "QS| GAPW|      max_l_rho0 :", &
               qs_control%gapw_control%lmax_rho0
            IF (qs_control%gapw_control%non_paw_atoms) THEN
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "QS| GAPW|      At least one kind is NOT PAW, i.e. it has only soft AO "
            END IF
            IF (qs_control%gapw_control%nopaw_as_gpw) THEN
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "QS| GAPW|      The NOT PAW atoms are treated fully GPW"
            END IF
         END IF
         IF (qs_control%gapw_xc) THEN
            SELECT CASE (qs_control%gapw_control%basis_1c)
            CASE (gapw_1c_orb)
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "QS| GAPW_XC|      One center basis from orbital basis primitives"
            CASE (gapw_1c_small)
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "QS| GAPW_XC|      One center basis extended with primitives (small:s)"
            CASE (gapw_1c_medium)
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "QS| GAPW_XC|      One center basis extended with primitives (medium:sp)"
            CASE (gapw_1c_large)
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "QS| GAPW_XC|      One center basis extended with primitives (large:spd)"
            CASE (gapw_1c_very_large)
               WRITE (unit=output_unit, fmt="(T2,A)") &
                  "QS| GAPW_XC|      One center basis extended with primitives (very large:spdf)"
            CASE DEFAULT
               cpabort("basis_1c incorrect")
            END SELECT
            WRITE (unit=output_unit, fmt="(T2,A,T73,ES8.1)") &
               "QS| GAPW_XC|                eps_fit:", &
               qs_control%gapw_control%eps_fit, &
               "QS| GAPW_XC|                eps_iso:", &
               qs_control%gapw_control%eps_iso, &
               "QS| GAPW_XC|                eps_svd:", &
               qs_control%gapw_control%eps_svd
            WRITE (unit=output_unit, fmt="(T2,A,T55,A30)") &
               "QS| GAPW_XC|atom-r-grid: quadrature:", &
               enum_i2c(enum, qs_control%gapw_control%quadrature)
            WRITE (unit=output_unit, fmt="(T2,A,T71,I10)") &
               "QS| GAPW_XC|   atom-s-grid:  max l :", &
               qs_control%gapw_control%lmax_sphere
         END IF
         IF (qs_control%mulliken_restraint) THEN
            WRITE (unit=output_unit, fmt="(T2,A,T73,ES8.1)") &
               "QS| Mulliken restraint target", qs_control%mulliken_restraint_control%target
            WRITE (unit=output_unit, fmt="(T2,A,T73,ES8.1)") &
               "QS| Mulliken restraint strength", qs_control%mulliken_restraint_control%strength
            WRITE (unit=output_unit, fmt="(T2,A,T73,I8)") &
               "QS| Mulliken restraint atoms: ", qs_control%mulliken_restraint_control%natoms
            WRITE (unit=output_unit, fmt="(5I8)") qs_control%mulliken_restraint_control%atoms
         END IF
         IF (qs_control%ddapc_restraint) THEN
            DO i = 1, SIZE(qs_control%ddapc_restraint_control)
               ddapc_restraint_control => qs_control%ddapc_restraint_control(i)
               IF (SIZE(qs_control%ddapc_restraint_control) > 1) &
                  WRITE (unit=output_unit, fmt="(T2,A,T3,I8)") &
                  "QS| parameters for DDAPC restraint number", i
               WRITE (unit=output_unit, fmt="(T2,A,T73,ES8.1)") &
                  "QS| ddapc restraint target", ddapc_restraint_control%target
               WRITE (unit=output_unit, fmt="(T2,A,T73,ES8.1)") &
                  "QS| ddapc restraint strength", ddapc_restraint_control%strength
               WRITE (unit=output_unit, fmt="(T2,A,T73,I8)") &
                  "QS| ddapc restraint atoms: ", ddapc_restraint_control%natoms
               WRITE (unit=output_unit, fmt="(5I8)") ddapc_restraint_control%atoms
               WRITE (unit=output_unit, fmt="(T2,A)") "Coefficients:"
               WRITE (unit=output_unit, fmt="(5F6.2)") ddapc_restraint_control%coeff
               SELECT CASE (ddapc_restraint_control%functional_form)
               CASE (do_ddapc_restraint)
                  WRITE (unit=output_unit, fmt="(T2,A,T61,A20)") &
                     "QS| ddapc restraint functional form :", "RESTRAINT"
               CASE (do_ddapc_constraint)
                  WRITE (unit=output_unit, fmt="(T2,A,T61,A20)") &
                     "QS| ddapc restraint functional form :", "CONSTRAINT"
               CASE DEFAULT
                  cpabort("Unknown ddapc restraint")
               END SELECT
            END DO
         END IF
         IF (qs_control%s2_restraint) THEN
            WRITE (unit=output_unit, fmt="(T2,A,T73,ES8.1)") &
               "QS| s2 restraint target", qs_control%s2_restraint_control%target
            WRITE (unit=output_unit, fmt="(T2,A,T73,ES8.1)") &
               "QS| s2 restraint strength", qs_control%s2_restraint_control%strength
            SELECT CASE (qs_control%s2_restraint_control%functional_form)
            CASE (do_s2_restraint)
               WRITE (unit=output_unit, fmt="(T2,A,T61,A20)") &
                  "QS| s2 restraint functional form :", "RESTRAINT"
               cpabort("Not yet implemented")
            CASE (do_s2_constraint)
               WRITE (unit=output_unit, fmt="(T2,A,T61,A20)") &
                  "QS| s2 restraint functional form :", "CONSTRAINT"
            CASE DEFAULT
               cpabort("Unknown ddapc restraint")
            END SELECT
         END IF
      END IF
      CALL cp_print_key_finished_output(output_unit, logger, print_section_vals, &
                                        "DFT_CONTROL_PARAMETERS")
 
      CALL timestop(handle)
 
   END SUBROUTINE write_qs_control
 
! **************************************************************************************************
!> \brief reads the input parameters needed for ddapc.
!> \param qs_control ...
!> \param qs_section ...
!> \param ddapc_restraint_section ...
!> \author fschiff
!> \note
!>      either reads DFT%QS%DDAPC_RESTRAINT or PROPERTIES%ET_coupling
!>      if(qs_section is present the DFT part is read, if ddapc_restraint_section
!>      is present ET_COUPLING is read. Avoid having both!!!
! **************************************************************************************************
   SUBROUTINE read_ddapc_section(qs_control, qs_section, ddapc_restraint_section)
 
      TYPE(qs_control_type), INTENT(INOUT)               :: qs_control
      TYPE(section_vals_type), OPTIONAL, POINTER         :: qs_section, ddapc_restraint_section
 
      INTEGER                                            :: i, j, jj, k, n_rep
      INTEGER, DIMENSION(:), POINTER                     :: tmplist
      REAL(kind=dp), DIMENSION(:), POINTER               :: rtmplist
      TYPE(ddapc_restraint_type), POINTER                :: ddapc_restraint_control
      TYPE(section_vals_type), POINTER                   :: ddapc_section
 
      IF (PRESENT(ddapc_restraint_section)) THEN
         IF (ASSOCIATED(qs_control%ddapc_restraint_control)) THEN
            IF (SIZE(qs_control%ddapc_restraint_control) >= 2) &
               cpabort("ET_COUPLING cannot be used in combination with a normal restraint")
         ELSE
            ddapc_section => ddapc_restraint_section
            ALLOCATE (qs_control%ddapc_restraint_control(1))
         END IF
      END IF
 
      IF (PRESENT(qs_section)) THEN
         NULLIFY (ddapc_section)
         ddapc_section => section_vals_get_subs_vals(qs_section, &
                                                     "DDAPC_RESTRAINT")
      END IF
 
      DO i = 1, SIZE(qs_control%ddapc_restraint_control)
 
         CALL ddapc_control_create(qs_control%ddapc_restraint_control(i))
         ddapc_restraint_control => qs_control%ddapc_restraint_control(i)
 
         CALL section_vals_val_get(ddapc_section, "STRENGTH", i_rep_section=i, &
                                   r_val=ddapc_restraint_control%strength)
         CALL section_vals_val_get(ddapc_section, "TARGET", i_rep_section=i, &
                                   r_val=ddapc_restraint_control%target)
         CALL section_vals_val_get(ddapc_section, "FUNCTIONAL_FORM", i_rep_section=i, &
                                   i_val=ddapc_restraint_control%functional_form)
         CALL section_vals_val_get(ddapc_section, "ATOMS", i_rep_section=i, &
                                   n_rep_val=n_rep)
         CALL section_vals_val_get(ddapc_section, "TYPE_OF_DENSITY", i_rep_section=i, &
                                   i_val=ddapc_restraint_control%density_type)
 
         jj = 0
         DO k = 1, n_rep
            CALL section_vals_val_get(ddapc_section, "ATOMS", i_rep_section=i, &
                                      i_rep_val=k, i_vals=tmplist)
            DO j = 1, SIZE(tmplist)
               jj = jj + 1
            END DO
         END DO
         IF (jj < 1) cpabort("Need at least 1 atom to use ddapc constraints")
         ddapc_restraint_control%natoms = jj
         IF (ASSOCIATED(ddapc_restraint_control%atoms)) &
            DEALLOCATE (ddapc_restraint_control%atoms)
         ALLOCATE (ddapc_restraint_control%atoms(ddapc_restraint_control%natoms))
         jj = 0
         DO k = 1, n_rep
            CALL section_vals_val_get(ddapc_section, "ATOMS", i_rep_section=i, &
                                      i_rep_val=k, i_vals=tmplist)
            DO j = 1, SIZE(tmplist)
               jj = jj + 1
               ddapc_restraint_control%atoms(jj) = tmplist(j)
            END DO
         END DO
 
         IF (ASSOCIATED(ddapc_restraint_control%coeff)) &
            DEALLOCATE (ddapc_restraint_control%coeff)
         ALLOCATE (ddapc_restraint_control%coeff(ddapc_restraint_control%natoms))
         ddapc_restraint_control%coeff = 1.0_dp
 
         CALL section_vals_val_get(ddapc_section, "COEFF", i_rep_section=i, &
                                   n_rep_val=n_rep)
         jj = 0
         DO k = 1, n_rep
            CALL section_vals_val_get(ddapc_section, "COEFF", i_rep_section=i, &
                                      i_rep_val=k, r_vals=rtmplist)
            DO j = 1, SIZE(rtmplist)
               jj = jj + 1
               IF (jj > ddapc_restraint_control%natoms) &
                  cpabort("Need the same number of coeff as there are atoms ")
               ddapc_restraint_control%coeff(jj) = rtmplist(j)
            END DO
         END DO
         IF (jj < ddapc_restraint_control%natoms .AND. jj /= 0) &
            cpabort("Need no or the same number of coeff as there are atoms.")
      END DO
      k = 0
      DO i = 1, SIZE(qs_control%ddapc_restraint_control)
         IF (qs_control%ddapc_restraint_control(i)%functional_form == &
             do_ddapc_constraint) k = k + 1
      END DO
      IF (k == 2) CALL cp_abort(__location__, &
                                "Only a single constraint possible yet, try to use restraints instead ")
 
   END SUBROUTINE read_ddapc_section
 
! **************************************************************************************************
!> \brief ...
!> \param dft_control ...
!> \param efield_section ...
!> \param cell ...
! **************************************************************************************************
   SUBROUTINE read_efield_sections(dft_control, efield_section, cell)
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(section_vals_type), POINTER                   :: efield_section
      TYPE(cell_type), OPTIONAL, POINTER                 :: cell
 
      CHARACTER(len=default_path_length)                 :: file_name
      INTEGER                                            :: i, io, j, n, unit_nr
      LOGICAL                                            :: amplitude_explicit, intensity_explicit
      REAL(kind=dp), DIMENSION(:), POINTER               :: tmp_vals
      TYPE(efield_type), POINTER                         :: efield
      TYPE(section_vals_type), POINTER                   :: tmp_section
 
      DO i = 1, SIZE(dft_control%efield_fields)
         NULLIFY (dft_control%efield_fields(i)%efield)
         ALLOCATE (dft_control%efield_fields(i)%efield)
         efield => dft_control%efield_fields(i)%efield
         NULLIFY (efield%envelop_i_vars, efield%envelop_r_vars)
         CALL section_vals_val_get(efield_section, "INTENSITY", i_rep_section=i, &
                                   r_val=efield%strength, explicit=intensity_explicit)
         CALL section_vals_val_get(efield_section, "AMPLITUDE", i_rep_section=i, &
                                   r_val=efield%amplitude, explicit=amplitude_explicit)
 
         IF (intensity_explicit .AND. amplitude_explicit) THEN
            cpabort("Both INTENSITY and AMPLITUDE provided in EFIELD section.")
         END IF
 
         IF (intensity_explicit) THEN
            efield%amplitude = sqrt(efield%strength/(3.50944_dp*10.0_dp**16))
         END IF
 
         IF (amplitude_explicit) THEN
            efield%strength = (3.50944_dp*10.0_dp**16)*(efield%amplitude)**2
         END IF
 
         CALL section_vals_val_get(efield_section, "POLARISATION", i_rep_section=i, &
                                   r_vals=tmp_vals)
         ALLOCATE (efield%polarisation(SIZE(tmp_vals)))
         efield%polarisation = tmp_vals
         IF (PRESENT(cell)) THEN
            IF (ASSOCIATED(cell)) CALL cell_transform_input_cartesian(cell, efield%polarisation(1:3))
         END IF
         CALL section_vals_val_get(efield_section, "PHASE", i_rep_section=i, &
                                   r_val=efield%phase_offset)
         CALL section_vals_val_get(efield_section, "ENVELOP", i_rep_section=i, &
                                   i_val=efield%envelop_id)
         CALL section_vals_val_get(efield_section, "WAVELENGTH", i_rep_section=i, &
                                   r_val=efield%wavelength)
         CALL section_vals_val_get(efield_section, "VEC_POT_INITIAL", i_rep_section=i, &
                                   r_vals=tmp_vals)
         efield%vec_pot_initial = tmp_vals
         IF (PRESENT(cell)) THEN
            IF (ASSOCIATED(cell)) CALL cell_transform_input_cartesian(cell, efield%vec_pot_initial(1:3))
         END IF
 
         IF (efield%envelop_id == constant_env) THEN
            ALLOCATE (efield%envelop_i_vars(2))
            tmp_section => section_vals_get_subs_vals(efield_section, "CONSTANT_ENV", i_rep_section=i)
            CALL section_vals_val_get(tmp_section, "START_STEP", &
                                      i_val=efield%envelop_i_vars(1))
            CALL section_vals_val_get(tmp_section, "END_STEP", &
                                      i_val=efield%envelop_i_vars(2))
         ELSE IF (efield%envelop_id == gaussian_env) THEN
            ALLOCATE (efield%envelop_r_vars(2))
            tmp_section => section_vals_get_subs_vals(efield_section, "GAUSSIAN_ENV", i_rep_section=i)
            CALL section_vals_val_get(tmp_section, "T0", &
                                      r_val=efield%envelop_r_vars(1))
            CALL section_vals_val_get(tmp_section, "SIGMA", &
                                      r_val=efield%envelop_r_vars(2))
         ELSE IF (efield%envelop_id == ramp_env) THEN
            ALLOCATE (efield%envelop_i_vars(4))
            tmp_section => section_vals_get_subs_vals(efield_section, "RAMP_ENV", i_rep_section=i)
            CALL section_vals_val_get(tmp_section, "START_STEP_IN", &
                                      i_val=efield%envelop_i_vars(1))
            CALL section_vals_val_get(tmp_section, "END_STEP_IN", &
                                      i_val=efield%envelop_i_vars(2))
            CALL section_vals_val_get(tmp_section, "START_STEP_OUT", &
                                      i_val=efield%envelop_i_vars(3))
            CALL section_vals_val_get(tmp_section, "END_STEP_OUT", &
                                      i_val=efield%envelop_i_vars(4))
         ELSE IF (efield%envelop_id == custom_env) THEN
            tmp_section => section_vals_get_subs_vals(efield_section, "CUSTOM_ENV", i_rep_section=i)
            CALL section_vals_val_get(tmp_section, "EFIELD_FILE_NAME", c_val=file_name)
            CALL open_file(file_name=trim(file_name), file_action="READ", file_status="OLD", unit_number=unit_nr)
            !Determine the number of lines in file
            n = 0
            DO WHILE (.true.)
               READ (unit_nr, *, iostat=io)
               IF (io /= 0) EXIT
               n = n + 1
            END DO
            rewind(unit_nr)
            ALLOCATE (efield%envelop_r_vars(n + 1))
            !Store the timestep of the list in the first entry of the r_vars
            CALL section_vals_val_get(tmp_section, "TIMESTEP", r_val=efield%envelop_r_vars(1))
            !Read the file
            DO j = 2, n + 1
               READ (unit_nr, *) efield%envelop_r_vars(j)
               efield%envelop_r_vars(j) = cp_unit_to_cp2k(efield%envelop_r_vars(j), "volt/m")
            END DO
            CALL close_file(unit_nr)
         END IF
      END DO
   END SUBROUTINE read_efield_sections
 
! **************************************************************************************************
!> \brief reads the input parameters needed real time propagation
!> \param dft_control ...
!> \param rtp_section ...
!> \author fschiff
! **************************************************************************************************
   SUBROUTINE read_rtp_section(dft_control, rtp_section)
 
      TYPE(dft_control_type), INTENT(INOUT)              :: dft_control
      TYPE(section_vals_type), POINTER                   :: rtp_section
 
      INTEGER                                            :: i, j, n_elems
      INTEGER, DIMENSION(:), POINTER                     :: tmp
      LOGICAL                                            :: is_present, local_moment_possible
      TYPE(section_vals_type), POINTER                   :: proj_mo_section, subsection
 
      ALLOCATE (dft_control%rtp_control)
      CALL section_vals_val_get(rtp_section, "MAX_ITER", &
                                i_val=dft_control%rtp_control%max_iter)
      CALL section_vals_val_get(rtp_section, "MAT_EXP", &
                                i_val=dft_control%rtp_control%mat_exp)
      CALL section_vals_val_get(rtp_section, "ASPC_ORDER", &
                                i_val=dft_control%rtp_control%aspc_order)
      CALL section_vals_val_get(rtp_section, "EXP_ACCURACY", &
                                r_val=dft_control%rtp_control%eps_exp)
      CALL section_vals_val_get(rtp_section, "RTBSE%_SECTION_PARAMETERS_", &
                                i_val=dft_control%rtp_control%rtp_method)
      CALL section_vals_val_get(rtp_section, "RTBSE%RTBSE_HAMILTONIAN", &
                                i_val=dft_control%rtp_control%rtbse_ham)
      CALL section_vals_val_get(rtp_section, "PROPAGATOR", &
                                i_val=dft_control%rtp_control%propagator)
      CALL section_vals_val_get(rtp_section, "EPS_ITER", &
                                r_val=dft_control%rtp_control%eps_ener)
      CALL section_vals_val_get(rtp_section, "INITIAL_WFN", &
                                i_val=dft_control%rtp_control%initial_wfn)
      CALL section_vals_val_get(rtp_section, "HFX_BALANCE_IN_CORE", &
                                l_val=dft_control%rtp_control%hfx_redistribute)
      CALL section_vals_val_get(rtp_section, "APPLY_WFN_MIX_INIT_RESTART", &
                                l_val=dft_control%rtp_control%apply_wfn_mix_init_restart)
      CALL section_vals_val_get(rtp_section, "APPLY_DELTA_PULSE", &
                                l_val=dft_control%rtp_control%apply_delta_pulse)
      CALL section_vals_val_get(rtp_section, "APPLY_DELTA_PULSE_MAG", &
                                l_val=dft_control%rtp_control%apply_delta_pulse_mag)
      CALL section_vals_val_get(rtp_section, "VELOCITY_GAUGE", &
                                l_val=dft_control%rtp_control%velocity_gauge)
      CALL section_vals_val_get(rtp_section, "VG_COM_NL", &
                                l_val=dft_control%rtp_control%nl_gauge_transform)
      CALL section_vals_val_get(rtp_section, "PERIODIC", &
                                l_val=dft_control%rtp_control%periodic)
      CALL section_vals_val_get(rtp_section, "DENSITY_PROPAGATION", &
                                l_val=dft_control%rtp_control%linear_scaling)
      CALL section_vals_val_get(rtp_section, "MCWEENY_MAX_ITER", &
                                i_val=dft_control%rtp_control%mcweeny_max_iter)
      CALL section_vals_val_get(rtp_section, "ACCURACY_REFINEMENT", &
                                i_val=dft_control%rtp_control%acc_ref)
      CALL section_vals_val_get(rtp_section, "MCWEENY_EPS", &
                                r_val=dft_control%rtp_control%mcweeny_eps)
      CALL section_vals_val_get(rtp_section, "DELTA_PULSE_SCALE", &
                                r_val=dft_control%rtp_control%delta_pulse_scale)
      CALL section_vals_val_get(rtp_section, "DELTA_PULSE_DIRECTION", &
                                i_vals=tmp)
      dft_control%rtp_control%delta_pulse_direction = tmp
      CALL section_vals_val_get(rtp_section, "SC_CHECK_START", &
                                i_val=dft_control%rtp_control%sc_check_start)
      proj_mo_section => section_vals_get_subs_vals(rtp_section, "PRINT%PROJECTION_MO")
      CALL section_vals_get(proj_mo_section, explicit=is_present)
      IF (is_present) THEN
         IF (dft_control%rtp_control%linear_scaling) &
            CALL cp_abort(__location__, &
                          "You have defined a time dependent projection of mos, but "// &
                          "only the density matrix is propagated (DENSITY_PROPAGATION "// &
                          ".TRUE.). Please either use MO-based real time DFT or do not "// &
                          "define any PRINT%PROJECTION_MO section")
         dft_control%rtp_control%is_proj_mo = .true.
      ELSE
         dft_control%rtp_control%is_proj_mo = .false.
      END IF
      ! Moment trace
      local_moment_possible = (dft_control%rtp_control%rtp_method == rtp_method_bse) .OR. &
                              ((.NOT. dft_control%rtp_control%periodic) .AND. dft_control%rtp_control%linear_scaling)
      ! TODO : Implement for other moment operators
      subsection => section_vals_get_subs_vals(rtp_section, "PRINT%MOMENTS")
      CALL section_vals_get(subsection, explicit=is_present)
      ! Trigger the flag
      dft_control%rtp_control%save_local_moments = &
         is_present .OR. dft_control%rtp_control%save_local_moments
      IF (is_present .AND. (.NOT. local_moment_possible)) THEN
         CALL cp_abort(__location__, "Moments trace printing only "// &
                       "implemented in non-periodic systems in linear scaling. "// &
                       "Please use DFT%PRINT%MOMENTS for other printing.")
      END IF
      CALL section_vals_val_get(rtp_section, "PRINT%MOMENTS%REFERENCE", &
                                i_val=dft_control%rtp_control%moment_trace_ref_type)
      CALL section_vals_val_get(rtp_section, "PRINT%MOMENTS%REFERENCE_POINT", &
                                r_vals=dft_control%rtp_control%moment_trace_user_ref_point)
      ! Moment Fourier transform
      subsection => section_vals_get_subs_vals(rtp_section, "PRINT%MOMENTS_FT")
      CALL section_vals_get(subsection, explicit=is_present)
      ! Trigger the flag
      dft_control%rtp_control%save_local_moments = &
         is_present .OR. dft_control%rtp_control%save_local_moments
      IF (is_present .AND. (.NOT. local_moment_possible)) THEN
         ! Not implemented
         CALL cp_abort(__location__, "Moments Fourier transform printing "// &
                       "implemented only for non-periodic systems in linear scaling.")
      END IF
      ! General FT settings
      CALL section_vals_val_get(rtp_section, "FT%DAMPING", &
                                r_val=dft_control%rtp_control%ft_damping)
      CALL section_vals_val_get(rtp_section, "FT%START_TIME", &
                                r_val=dft_control%rtp_control%ft_t0)
      ! Padé settings
      subsection => section_vals_get_subs_vals(rtp_section, "FT%PADE")
      CALL section_vals_val_get(subsection, "_SECTION_PARAMETERS_", &
                                l_val=dft_control%rtp_control%pade_requested)
      CALL section_vals_val_get(subsection, "E_MIN", &
                                r_val=dft_control%rtp_control%pade_e_min)
      CALL section_vals_val_get(subsection, "E_STEP", &
                                r_val=dft_control%rtp_control%pade_e_step)
      CALL section_vals_val_get(subsection, "E_MAX", &
                                r_val=dft_control%rtp_control%pade_e_max)
      CALL section_vals_val_get(subsection, "FIT_E_MIN", &
                                r_val=dft_control%rtp_control%pade_fit_e_min)
      CALL section_vals_val_get(subsection, "FIT_E_MAX", &
                                r_val=dft_control%rtp_control%pade_fit_e_max)
      ! If default settings used for fit_e_min/max, rewrite with appropriate values
      IF (dft_control%rtp_control%pade_fit_e_min < 0) THEN
         dft_control%rtp_control%pade_fit_e_min = dft_control%rtp_control%pade_e_min
      END IF
      IF (dft_control%rtp_control%pade_fit_e_max < 0) THEN
         dft_control%rtp_control%pade_fit_e_max = dft_control%rtp_control%pade_e_max
      END IF
      ! Polarizability settings
      subsection => section_vals_get_subs_vals(rtp_section, "PRINT%POLARIZABILITY")
      CALL section_vals_get(subsection, explicit=is_present)
      ! Trigger the flag
      dft_control%rtp_control%save_local_moments = &
         is_present .OR. dft_control%rtp_control%save_local_moments
      IF (is_present .AND. (.NOT. local_moment_possible)) THEN
         ! Not implemented
         CALL cp_abort(__location__, "Polarizability printing "// &
                       "implemented only for non-periodic systems.")
      END IF
      CALL section_vals_val_get(subsection, "ELEMENT", explicit=is_present, n_rep_val=n_elems)
      NULLIFY (dft_control%rtp_control%print_pol_elements)
      IF (is_present) THEN
         ! Explicit list of elements
         ! Allocate the array
         ALLOCATE (dft_control%rtp_control%print_pol_elements(n_elems, 2))
         DO i = 1, n_elems
            CALL section_vals_val_get(subsection, "ELEMENT", i_vals=tmp, i_rep_val=i)
            dft_control%rtp_control%print_pol_elements(i, :) = tmp(:)
         END DO
         ! Do basic sanity checks for pol_element
         DO i = 1, n_elems
            DO j = 1, 2
               IF (dft_control%rtp_control%print_pol_elements(i, j) > 3 .OR. &
                   dft_control%rtp_control%print_pol_elements(i, j) < 1) &
                  cpabort("Polarisation tensor element not 1,2 or 3 in at least one index")
            END DO
         END DO
      END IF
 
      ! Finally, allow printing of FT observables also in the case when they are not explicitly
      ! required, but they are available, i.e. non-periodic linear scaling calculation
      dft_control%rtp_control%save_local_moments = &
         dft_control%rtp_control%save_local_moments .OR. &
         ((.NOT. dft_control%rtp_control%periodic) .AND. dft_control%rtp_control%linear_scaling)
 
   END SUBROUTINE read_rtp_section
! **************************************************************************************************
!> \brief Tries to guess the elements of polarization to print
!> \param dftc DFT parameters
!> \param elems 2D array, where the guessed element indeces are stored
!> \date 11.2025
!> \author Stepan Marek
! **************************************************************************************************
   SUBROUTINE guess_pol_elements(dftc, elems)
      TYPE(dft_control_type)                             :: dftc
      INTEGER, DIMENSION(:, :), POINTER                  :: elems
 
      INTEGER                                            :: i, i_nonzero, n_nonzero
      LOGICAL                                            :: pol_vector_known
      REAL(kind=dp), DIMENSION(3)                        :: pol_vector
 
      pol_vector_known = .false.
 
      ! TODO : More relevant elements for magnetic pulse?
      IF (dftc%rtp_control%apply_delta_pulse .OR. dftc%rtp_control%apply_delta_pulse_mag) THEN
         pol_vector(:) = real(dftc%rtp_control%delta_pulse_direction(:), kind=dp)
      ELSE
         ! Maybe RT field is applied?
         pol_vector(:) = dftc%efield_fields(1)%efield%polarisation(:)
      END IF
      IF (dot_product(pol_vector, pol_vector) > 0.0_dp) pol_vector_known = .true.
 
      IF (.NOT. pol_vector_known) THEN
         cpabort("Cannot guess polarization elements - please specify!")
      ELSE
         ! Check whether just one element is non-zero
         n_nonzero = 0
         DO i = 1, 3
            IF (pol_vector(i) /= 0.0_dp) THEN
               n_nonzero = n_nonzero + 1
               i_nonzero = i
            END IF
         END DO
         IF (n_nonzero > 1) THEN
            CALL cp_abort(__location__, &
                          "More than one non-zero field elements - "// &
                          "cannot guess polarizability elements - please specify!")
         ELSE IF (n_nonzero == 0) THEN
            CALL cp_abort(__location__, &
                          "No non-zero field elements - "// &
                          "cannot guess polarizability elements - please specify!")
         ELSE
            ! Clear guess can be made
            NULLIFY (elems)
            ALLOCATE (elems(3, 2))
            DO i = 1, 3
               elems(i, 1) = i
               elems(i, 2) = i_nonzero
            END DO
         END IF
      END IF
   END SUBROUTINE guess_pol_elements
 
! **************************************************************************************************
!> \brief Parses the BLOCK_LIST keywords from the ADMM section
!> \param admm_control ...
!> \param dft_section ...
! **************************************************************************************************
   SUBROUTINE read_admm_block_list(admm_control, dft_section)
      TYPE(admm_control_type), POINTER                   :: admm_control
      TYPE(section_vals_type), POINTER                   :: dft_section
 
      INTEGER                                            :: irep, list_size, n_rep
      INTEGER, DIMENSION(:), POINTER                     :: tmplist
 
      NULLIFY (tmplist)
 
      CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%BLOCK_LIST", &
                                n_rep_val=n_rep)
 
      ALLOCATE (admm_control%blocks(n_rep))
 
      DO irep = 1, n_rep
         CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%BLOCK_LIST", &
                                   i_rep_val=irep, i_vals=tmplist)
         list_size = SIZE(tmplist)
         ALLOCATE (admm_control%blocks(irep)%list(list_size))
         admm_control%blocks(irep)%list(:) = tmplist(:)
      END DO
 
   END SUBROUTINE read_admm_block_list
 
! **************************************************************************************************
!> \brief ...
!> \param dft_control ...
!> \param hairy_probes_section ...
!> \param
!> \param
! **************************************************************************************************
   SUBROUTINE read_hairy_probes_sections(dft_control, hairy_probes_section)
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(section_vals_type), POINTER                   :: hairy_probes_section
 
      INTEGER                                            :: i, j, jj, kk, n_rep
      INTEGER, DIMENSION(:), POINTER                     :: tmplist
 
      DO i = 1, SIZE(dft_control%probe)
         NULLIFY (dft_control%probe(i)%atom_ids)
 
         CALL section_vals_val_get(hairy_probes_section, "ATOM_IDS", i_rep_section=i, n_rep_val=n_rep)
         jj = 0
         DO kk = 1, n_rep
            CALL section_vals_val_get(hairy_probes_section, "ATOM_IDS", i_rep_section=i, i_rep_val=kk, i_vals=tmplist)
            jj = jj + SIZE(tmplist)
         END DO
 
         dft_control%probe(i)%natoms = jj
         IF (dft_control%probe(i)%natoms < 1) &
            cpabort("Need at least 1 atom to use hair probes formalism")
         ALLOCATE (dft_control%probe(i)%atom_ids(dft_control%probe(i)%natoms))
 
         jj = 0
         DO kk = 1, n_rep
            CALL section_vals_val_get(hairy_probes_section, "ATOM_IDS", i_rep_section=i, i_rep_val=kk, i_vals=tmplist)
            DO j = 1, SIZE(tmplist)
               jj = jj + 1
               dft_control%probe(i)%atom_ids(jj) = tmplist(j)
            END DO
         END DO
 
         CALL section_vals_val_get(hairy_probes_section, "MU", i_rep_section=i, r_val=dft_control%probe(i)%mu)
 
         CALL section_vals_val_get(hairy_probes_section, "T", i_rep_section=i, r_val=dft_control%probe(i)%T)
 
         CALL section_vals_val_get(hairy_probes_section, "ALPHA", i_rep_section=i, r_val=dft_control%probe(i)%alpha)
 
         CALL section_vals_val_get(hairy_probes_section, "eps_hp", i_rep_section=i, r_val=dft_control%probe(i)%eps_hp)
      END DO
 
   END SUBROUTINE read_hairy_probes_sections
! **************************************************************************************************
 
END MODULE cp_control_utils