LCOV - code coverage report
Current view: top level - src - cp_control_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:24d69ee) Lines: 86.2 % 1916 1651
Test Date: 2026-09-03 07:32:15 Functions: 100.0 % 17 17

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief Utilities to set up the control types
      10              : ! **************************************************************************************************
      11              : MODULE cp_control_utils
      12              :    USE bibliography,                    ONLY: &
      13              :         Andreussi2012, Andreussi2019, Chai2025a, Dewar1977, Dewar1985, Elstner1998, Fattebert2002, &
      14              :         Grimme2017, Hu2007, Katbashev2025, Krack2000, Lippert1997, Lippert1999, Porezag1995, &
      15              :         Pracht2019, Repasky2002, Rocha2006, Schenter2008, Seifert1996, Souza2002, Stengel2009, &
      16              :         Stewart1989, Stewart2007, Thiel1992, Umari2002, VanVoorhis2015, VandeVondele2005a, &
      17              :         VandeVondele2005b, Yin2017, Zhechkov2005, cite_reference
      18              :    USE cell_types,                      ONLY: cell_transform_input_cartesian,&
      19              :                                               cell_type
      20              :    USE cp_control_types,                ONLY: &
      21              :         admm_control_create, admm_control_type, ddapc_control_create, ddapc_restraint_type, &
      22              :         dft_control_create, dft_control_type, efield_type, expot_control_create, &
      23              :         maxwell_control_create, qs_control_type, rixs_control_type, tddfpt2_control_type, &
      24              :         xtb_control_type, xtb_reference_cli_type
      25              :    USE cp_files,                        ONLY: close_file,&
      26              :                                               open_file
      27              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      28              :                                               cp_logger_type
      29              :    USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
      30              :                                               cp_print_key_unit_nr
      31              :    USE cp_parser_methods,               ONLY: parser_read_line
      32              :    USE cp_parser_types,                 ONLY: cp_parser_type,&
      33              :                                               parser_create,&
      34              :                                               parser_release,&
      35              :                                               parser_reset
      36              :    USE cp_spline_utils,                 ONLY: pw_interp
      37              :    USE cp_units,                        ONLY: cp_unit_from_cp2k,&
      38              :                                               cp_unit_to_cp2k
      39              :    USE eeq_input,                       ONLY: read_eeq_param
      40              :    USE force_fields_input,              ONLY: read_gp_section
      41              :    USE input_constants,                 ONLY: &
      42              :         admm1_type, admm2_type, admmp_type, admmq_type, admms_type, constant_env, custom_env, &
      43              :         do_admm_aux_exch_func_bee, do_admm_aux_exch_func_bee_libxc, do_admm_aux_exch_func_default, &
      44              :         do_admm_aux_exch_func_default_libxc, do_admm_aux_exch_func_none, &
      45              :         do_admm_aux_exch_func_opt, do_admm_aux_exch_func_opt_libxc, do_admm_aux_exch_func_pbex, &
      46              :         do_admm_aux_exch_func_pbex_libxc, do_admm_aux_exch_func_sx_libxc, &
      47              :         do_admm_basis_projection, do_admm_blocked_projection, do_admm_blocking_purify_full, &
      48              :         do_admm_charge_constrained_projection, do_admm_exch_scaling_merlot, &
      49              :         do_admm_exch_scaling_none, do_admm_purify_cauchy, do_admm_purify_cauchy_subspace, &
      50              :         do_admm_purify_mcweeny, do_admm_purify_mo_diag, do_admm_purify_mo_no_diag, &
      51              :         do_admm_purify_none, do_admm_purify_none_dm, do_ddapc_constraint, do_ddapc_restraint, &
      52              :         do_method_am1, do_method_dftb, do_method_gapw, do_method_gapw_xc, do_method_gpw, &
      53              :         do_method_lrigpw, do_method_mndo, do_method_mndod, do_method_ofgpw, do_method_pdg, &
      54              :         do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rigpw, &
      55              :         do_method_rm1, do_method_xtb, do_pwgrid_ns_fullspace, do_pwgrid_ns_halfspace, &
      56              :         do_pwgrid_spherical, do_s2_constraint, do_s2_restraint, do_se_is_kdso, do_se_is_kdso_d, &
      57              :         do_se_is_slater, do_se_lr_ewald, do_se_lr_ewald_gks, do_se_lr_ewald_r3, do_se_lr_none, &
      58              :         gapw_1c_large, gapw_1c_medium, gapw_1c_orb, gapw_1c_small, gapw_1c_very_large, &
      59              :         gaussian_env, gfn1xtb, gfn_tblite, kg_tnadd_embed, kg_tnadd_embed_ri, no_admm_type, &
      60              :         numerical, ramp_env, real_time_propagation, rtp_method_bse, rtp_method_bse_linearized, &
      61              :         sccs_andreussi, sccs_derivative_cd3, sccs_derivative_cd5, sccs_derivative_cd7, &
      62              :         sccs_derivative_fft, sccs_fattebert_gygi, sccs_saa_andreussi, sic_ad, sic_eo, &
      63              :         sic_list_all, sic_list_unpaired, sic_mauri_spz, sic_mauri_us, sic_none, slater, &
      64              :         tblite_cli_born_kernel_auto, tblite_cli_solution_state_gsolv, tblite_cli_solvation_alpb, &
      65              :         tblite_cli_solvation_cpcm, tblite_cli_solvation_gb, tblite_cli_solvation_gbe, &
      66              :         tblite_cli_solvation_gbsa, tblite_guess_ceh, tblite_mixer_memory_inherit, &
      67              :         tblite_scc_mixer_auto, tblite_scc_mixer_cp2k, tblite_scc_mixer_none, &
      68              :         tblite_scc_mixer_tblite, tblite_solver_gvd, tblite_solver_gvr, tddfpt_dipole_length, &
      69              :         tddfpt_kernel_stda, use_mom_ref_user, xtb_vdw_type_d3, xtb_vdw_type_d4, xtb_vdw_type_none
      70              :    USE input_cp2k_check,                ONLY: xc_functionals_expand
      71              :    USE input_cp2k_dft,                  ONLY: create_dft_section
      72              :    USE input_enumeration_types,         ONLY: enum_i2c,&
      73              :                                               enumeration_type
      74              :    USE input_keyword_types,             ONLY: keyword_get,&
      75              :                                               keyword_type
      76              :    USE input_section_types,             ONLY: &
      77              :         section_get_ival, section_get_keyword, section_release, section_type, section_vals_get, &
      78              :         section_vals_get_subs_vals, section_vals_type, section_vals_val_get, section_vals_val_set
      79              :    USE kinds,                           ONLY: default_path_length,&
      80              :                                               default_string_length,&
      81              :                                               dp
      82              :    USE mathconstants,                   ONLY: fourpi
      83              :    USE pair_potential_types,            ONLY: pair_potential_reallocate
      84              :    USE periodic_table,                  ONLY: get_ptable_info
      85              :    USE qs_cdft_utils,                   ONLY: read_cdft_control_section
      86              :    USE smeagol_control_types,           ONLY: read_smeagol_control
      87              :    USE string_utilities,                ONLY: uppercase
      88              :    USE util,                            ONLY: sort
      89              :    USE xas_tdp_types,                   ONLY: read_xas_tdp_control
      90              :    USE xc,                              ONLY: xc_uses_kinetic_energy_density,&
      91              :                                               xc_uses_norm_drho
      92              :    USE xc_input_constants,              ONLY: xc_deriv_collocate
      93              :    USE xc_write_output,                 ONLY: xc_write
      94              : #include "./base/base_uses.f90"
      95              : 
      96              :    IMPLICIT NONE
      97              : 
      98              :    PRIVATE
      99              : 
     100              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_control_utils'
     101              : 
     102              :    PUBLIC :: read_dft_control, &
     103              :              read_rixs_control, &
     104              :              read_mgrid_section, &
     105              :              read_qs_section, &
     106              :              read_tddfpt2_control, &
     107              :              write_dft_control, &
     108              :              write_qs_control, &
     109              :              write_admm_control, &
     110              :              read_ddapc_section
     111              : CONTAINS
     112              : 
     113              : ! **************************************************************************************************
     114              : !> \brief ...
     115              : !> \param dft_control ...
     116              : !> \param dft_section ...
     117              : !> \param cell ...
     118              : ! **************************************************************************************************
     119       162576 :    SUBROUTINE read_dft_control(dft_control, dft_section, cell)
     120              :       TYPE(dft_control_type), POINTER                    :: dft_control
     121              :       TYPE(section_vals_type), POINTER                   :: dft_section
     122              :       TYPE(cell_type), OPTIONAL, POINTER                 :: cell
     123              : 
     124              :       CHARACTER(len=default_path_length)                 :: basis_set_file_name, gauxc_model_name, &
     125              :                                                             intensities_file_name, &
     126              :                                                             potential_file_name
     127              :       CHARACTER(LEN=default_string_length), &
     128         9032 :          DIMENSION(:), POINTER                           :: tmpstringlist
     129              :       INTEGER                                            :: admmtype, irep, isize, kg_tnadd_method, &
     130              :                                                             method_id, nrep, xc_deriv_method_id
     131              :       LOGICAL :: at_end, do_hfx, do_ot, do_rpa_admm, do_rtp, exopt1, exopt2, exopt3, explicit, &
     132              :          is_present, l_param, local_moment_possible, native_skala_grid, not_SE, was_present
     133              :       REAL(KIND=dp)                                      :: density_cut, gradient_cut, tau_cut
     134         9032 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: pol
     135              :       TYPE(cp_logger_type), POINTER                      :: logger
     136              :       TYPE(cp_parser_type)                               :: parser
     137              :       TYPE(section_vals_type), POINTER :: hairy_probes_section, hfx_section, kg_xc_fun_section, &
     138              :          kg_xc_section, maxwell_section, sccs_section, scf_section, tmp_section, xc_fun_section, &
     139              :          xc_gauxc_subsection, xc_section
     140              : 
     141         9032 :       was_present = .FALSE.
     142              : 
     143         9032 :       logger => cp_get_default_logger()
     144              : 
     145         9032 :       NULLIFY (kg_xc_fun_section, kg_xc_section, tmp_section, xc_fun_section, xc_section, xc_gauxc_subsection)
     146        45160 :       ALLOCATE (dft_control)
     147         9032 :       CALL dft_control_create(dft_control)
     148              :       ! determine wheather this is a semiempirical or DFTB run
     149              :       ! --> (no XC section needs to be provided)
     150         9032 :       not_SE = .TRUE.
     151         9032 :       CALL section_vals_val_get(dft_section, "QS%METHOD", i_val=method_id)
     152         2526 :       SELECT CASE (method_id)
     153              :       CASE (do_method_dftb, do_method_xtb, do_method_mndo, do_method_am1, do_method_pm3, do_method_pnnl, &
     154              :             do_method_pm6, do_method_pm6fm, do_method_pdg, do_method_rm1, do_method_mndod)
     155         9032 :          not_SE = .FALSE.
     156              :       END SELECT
     157              :       ! Check for XC section and XC_FUNCTIONAL section
     158         9032 :       xc_section => section_vals_get_subs_vals(dft_section, "XC")
     159         9032 :       CALL section_vals_get(xc_section, explicit=is_present)
     160         9032 :       IF (.NOT. is_present .AND. not_SE) THEN
     161            0 :          CPABORT("XC section missing.")
     162              :       END IF
     163         9032 :       IF (is_present) THEN
     164         6522 :          CALL section_vals_val_get(xc_section, "density_cutoff", r_val=density_cut)
     165         6522 :          CALL section_vals_val_get(xc_section, "gradient_cutoff", r_val=gradient_cut)
     166         6522 :          CALL section_vals_val_get(xc_section, "tau_cutoff", r_val=tau_cut)
     167              :          ! Perform numerical stability checks and possibly correct the issues
     168         6522 :          IF (density_cut <= EPSILON(0.0_dp)*100.0_dp) THEN
     169              :             CALL cp_warn(__LOCATION__, &
     170              :                          "DENSITY_CUTOFF lower than 100*EPSILON, where EPSILON is the machine precision. "// &
     171            0 :                          "This may lead to numerical problems. Setting up shake_tol to 100*EPSILON! ")
     172              :          END IF
     173         6522 :          density_cut = MAX(EPSILON(0.0_dp)*100.0_dp, density_cut)
     174         6522 :          IF (gradient_cut <= EPSILON(0.0_dp)*100.0_dp) THEN
     175              :             CALL cp_warn(__LOCATION__, &
     176              :                          "GRADIENT_CUTOFF lower than 100*EPSILON, where EPSILON is the machine precision. "// &
     177            0 :                          "This may lead to numerical problems. Setting up shake_tol to 100*EPSILON! ")
     178              :          END IF
     179         6522 :          gradient_cut = MAX(EPSILON(0.0_dp)*100.0_dp, gradient_cut)
     180         6522 :          IF (tau_cut <= EPSILON(0.0_dp)*100.0_dp) THEN
     181              :             CALL cp_warn(__LOCATION__, &
     182              :                          "TAU_CUTOFF lower than 100*EPSILON, where EPSILON is the machine precision. "// &
     183            0 :                          "This may lead to numerical problems. Setting up shake_tol to 100*EPSILON! ")
     184              :          END IF
     185         6522 :          tau_cut = MAX(EPSILON(0.0_dp)*100.0_dp, tau_cut)
     186         6522 :          CALL section_vals_val_set(xc_section, "density_cutoff", r_val=density_cut)
     187         6522 :          CALL section_vals_val_set(xc_section, "gradient_cutoff", r_val=gradient_cut)
     188         6522 :          CALL section_vals_val_set(xc_section, "tau_cutoff", r_val=tau_cut)
     189              :       END IF
     190         9032 :       xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
     191         9032 :       CALL section_vals_get(xc_fun_section, explicit=is_present)
     192         9032 :       IF (.NOT. is_present .AND. not_SE) THEN
     193            0 :          CPABORT("XC_FUNCTIONAL section missing.")
     194              :       END IF
     195              : 
     196         9032 :       dft_control%use_gauxc = .FALSE.
     197         9032 :       IF (is_present) THEN
     198         6522 :          xc_gauxc_subsection => section_vals_get_subs_vals(xc_fun_section, "GAUXC")
     199         6522 :          CALL section_vals_get(xc_gauxc_subsection, explicit=dft_control%use_gauxc)
     200              :       END IF
     201              : 
     202         9032 :       scf_section => section_vals_get_subs_vals(dft_section, "SCF")
     203         9032 :       CALL section_vals_val_get(dft_section, "UKS", l_val=dft_control%uks)
     204         9032 :       CALL section_vals_val_get(dft_section, "ROKS", l_val=dft_control%roks)
     205         9032 :       IF (dft_control%uks .OR. dft_control%roks) THEN
     206         1903 :          dft_control%nspins = 2
     207              :       ELSE
     208         7129 :          dft_control%nspins = 1
     209              :       END IF
     210              : 
     211         9032 :       dft_control%lsd = (dft_control%nspins > 1)
     212         9032 :       dft_control%use_kinetic_energy_density = xc_uses_kinetic_energy_density(xc_fun_section, dft_control%lsd)
     213         9032 :       IF (dft_control%use_gauxc) THEN
     214              :          native_skala_grid = .FALSE.
     215          142 :          CALL section_vals_val_get(xc_gauxc_subsection, "NATIVE_GRID", l_val=native_skala_grid)
     216          142 :          CALL section_vals_val_get(xc_gauxc_subsection, "MODEL", c_val=gauxc_model_name)
     217          142 :          gauxc_model_name = ADJUSTL(gauxc_model_name)
     218          142 :          CALL uppercase(gauxc_model_name)
     219          142 :          IF (native_skala_grid .OR. &
     220              :              (TRIM(gauxc_model_name) /= "" .AND. TRIM(gauxc_model_name) /= "NONE")) THEN
     221          134 :             dft_control%use_kinetic_energy_density = .TRUE.
     222              :          END IF
     223              :       END IF
     224         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "KG_METHOD")
     225         9032 :       CALL section_vals_get(tmp_section, explicit=explicit)
     226         9032 :       IF (explicit) THEN
     227           82 :          CALL section_vals_val_get(tmp_section, "TNADD_METHOD", i_val=kg_tnadd_method)
     228           82 :          IF (kg_tnadd_method == kg_tnadd_embed .OR. kg_tnadd_method == kg_tnadd_embed_ri) THEN
     229           64 :             kg_xc_section => section_vals_get_subs_vals(tmp_section, "XC")
     230           64 :             kg_xc_fun_section => section_vals_get_subs_vals(kg_xc_section, "XC_FUNCTIONAL")
     231           64 :             CALL section_vals_get(kg_xc_fun_section, explicit=is_present)
     232           64 :             IF (is_present) THEN
     233              :                dft_control%use_kinetic_energy_density = dft_control%use_kinetic_energy_density .OR. &
     234              :                                                         xc_uses_kinetic_energy_density(kg_xc_fun_section, &
     235           68 :                                                                                        dft_control%lsd)
     236              :             END IF
     237              :          END IF
     238              :       END IF
     239              : 
     240         9032 :       xc_deriv_method_id = section_get_ival(xc_section, "XC_GRID%XC_DERIV")
     241              :       dft_control%drho_by_collocation = (xc_uses_norm_drho(xc_fun_section, dft_control%lsd) &
     242         9032 :                                          .AND. (xc_deriv_method_id == xc_deriv_collocate))
     243         9032 :       IF (dft_control%drho_by_collocation) THEN
     244            0 :          CPABORT("derivatives by collocation not implemented")
     245              :       END IF
     246              : 
     247              :       ! Automatic auxiliary basis set generation
     248         9032 :       CALL section_vals_val_get(dft_section, "AUTO_BASIS", n_rep_val=nrep)
     249        18064 :       DO irep = 1, nrep
     250         9032 :          CALL section_vals_val_get(dft_section, "AUTO_BASIS", i_rep_val=irep, c_vals=tmpstringlist)
     251        18064 :          IF (SIZE(tmpstringlist) == 2) THEN
     252         9032 :             CALL uppercase(tmpstringlist(2))
     253        17908 :             SELECT CASE (tmpstringlist(2))
     254              :             CASE ("X")
     255         8876 :                SELECT CASE (tmpstringlist(1))
     256              :                CASE ("X")
     257              :                   ! Do nothing
     258              :                CASE DEFAULT
     259              :                   CALL cp_abort(__LOCATION__, &
     260              :                                 "AUTO_BASIS: the size <X> is invalid for the "// &
     261              :                                 "type <"//TRIM(ADJUSTL(tmpstringlist(1)))//">; "// &
     262              :                                 "use one of SMALL, MEDIUM, LARGE, HUGE for "// &
     263              :                                 "the size. The syntax AUTO_BASIS X X is a "// &
     264              :                                 "reserved case for using NO automatically "// &
     265         8876 :                                 "generated basis sets.")
     266              :                END SELECT
     267              :             CASE ("SMALL")
     268           54 :                isize = 0
     269              :             CASE ("MEDIUM")
     270           54 :                isize = 1
     271              :             CASE ("LARGE")
     272            0 :                isize = 2
     273              :             CASE ("HUGE")
     274            8 :                isize = 3
     275              :             CASE DEFAULT
     276         9032 :                CPWARN("Unknown basis size in AUTO_BASIS keyword:"//TRIM(tmpstringlist(1)))
     277              :             END SELECT
     278              :             !
     279         9034 :             SELECT CASE (tmpstringlist(1))
     280              :             CASE ("X")
     281              :             CASE ("RI_AUX")
     282            2 :                dft_control%auto_basis_ri_aux = isize
     283              :             CASE ("AUX_FIT")
     284            0 :                dft_control%auto_basis_aux_fit = isize
     285              :             CASE ("LRI_AUX")
     286            0 :                dft_control%auto_basis_lri_aux = isize
     287              :             CASE ("P_LRI_AUX")
     288            0 :                dft_control%auto_basis_p_lri_aux = isize
     289              :             CASE ("RI_HXC")
     290            0 :                dft_control%auto_basis_ri_hxc = isize
     291              :             CASE ("RI_XAS")
     292           64 :                dft_control%auto_basis_ri_xas = isize
     293              :             CASE ("RI_HFX")
     294           90 :                dft_control%auto_basis_ri_hfx = isize
     295              :             CASE DEFAULT
     296         9032 :                CPWARN("Unknown basis type in AUTO_BASIS keyword:"//TRIM(tmpstringlist(1)))
     297              :             END SELECT
     298              :          ELSE
     299              :             CALL cp_abort(__LOCATION__, &
     300            0 :                           "AUTO_BASIS keyword in &DFT section has a wrong number of arguments.")
     301              :          END IF
     302              :       END DO
     303              : 
     304              :       !! check if we do wavefunction fitting
     305         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD")
     306         9032 :       CALL section_vals_get(tmp_section, explicit=is_present)
     307              :       !
     308         9032 :       hfx_section => section_vals_get_subs_vals(xc_section, "HF")
     309         9032 :       CALL section_vals_get(hfx_section, explicit=do_hfx)
     310         9032 :       CALL section_vals_val_get(xc_section, "WF_CORRELATION%RI_RPA%ADMM", l_val=do_rpa_admm)
     311         9032 :       is_present = is_present .AND. (do_hfx .OR. do_rpa_admm)
     312              :       !
     313         9032 :       dft_control%do_admm = is_present
     314         9032 :       dft_control%do_admm_mo = .FALSE.
     315         9032 :       dft_control%do_admm_dm = .FALSE.
     316         9032 :       IF (is_present) THEN
     317              :          do_ot = .FALSE.
     318          524 :          CALL section_vals_val_get(scf_section, "OT%_SECTION_PARAMETERS_", l_val=do_ot)
     319          524 :          CALL admm_control_create(dft_control%admm_control)
     320              : 
     321          524 :          CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%ADMM_TYPE", i_val=admmtype)
     322          524 :          CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%ADMM_PURIFICATION_METHOD", explicit=exopt1)
     323          524 :          CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%METHOD", explicit=exopt2)
     324          524 :          CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%EXCH_SCALING_MODEL", explicit=exopt3)
     325          524 :          dft_control%admm_control%admm_type = admmtype
     326          506 :          SELECT CASE (admmtype)
     327              :          CASE (no_admm_type)
     328          506 :             CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%ADMM_PURIFICATION_METHOD", i_val=method_id)
     329          506 :             dft_control%admm_control%purification_method = method_id
     330          506 :             CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%METHOD", i_val=method_id)
     331          506 :             dft_control%admm_control%method = method_id
     332          506 :             CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%EXCH_SCALING_MODEL", i_val=method_id)
     333          506 :             dft_control%admm_control%scaling_model = method_id
     334              :          CASE (admm1_type)
     335              :             ! METHOD BASIS_PROJECTION
     336              :             ! ADMM_PURIFICATION_METHOD choose
     337              :             ! EXCH_SCALING_MODEL NONE
     338            4 :             CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%ADMM_PURIFICATION_METHOD", i_val=method_id)
     339            4 :             dft_control%admm_control%purification_method = method_id
     340            4 :             dft_control%admm_control%method = do_admm_basis_projection
     341            4 :             dft_control%admm_control%scaling_model = do_admm_exch_scaling_none
     342              :          CASE (admm2_type)
     343              :             ! METHOD BASIS_PROJECTION
     344              :             ! ADMM_PURIFICATION_METHOD NONE
     345              :             ! EXCH_SCALING_MODEL NONE
     346            2 :             dft_control%admm_control%purification_method = do_admm_purify_none
     347            2 :             dft_control%admm_control%method = do_admm_basis_projection
     348            2 :             dft_control%admm_control%scaling_model = do_admm_exch_scaling_none
     349              :          CASE (admms_type)
     350              :             ! ADMM_PURIFICATION_METHOD NONE
     351              :             ! METHOD CHARGE_CONSTRAINED_PROJECTION
     352              :             ! EXCH_SCALING_MODEL MERLOT
     353            8 :             dft_control%admm_control%purification_method = do_admm_purify_none
     354            8 :             dft_control%admm_control%method = do_admm_charge_constrained_projection
     355            8 :             dft_control%admm_control%scaling_model = do_admm_exch_scaling_merlot
     356              :          CASE (admmp_type)
     357              :             ! ADMM_PURIFICATION_METHOD NONE
     358              :             ! METHOD BASIS_PROJECTION
     359              :             ! EXCH_SCALING_MODEL MERLOT
     360            2 :             dft_control%admm_control%purification_method = do_admm_purify_none
     361            2 :             dft_control%admm_control%method = do_admm_basis_projection
     362            2 :             dft_control%admm_control%scaling_model = do_admm_exch_scaling_merlot
     363              :          CASE (admmq_type)
     364              :             ! ADMM_PURIFICATION_METHOD NONE
     365              :             ! METHOD CHARGE_CONSTRAINED_PROJECTION
     366              :             ! EXCH_SCALING_MODEL NONE
     367            2 :             dft_control%admm_control%purification_method = do_admm_purify_none
     368            2 :             dft_control%admm_control%method = do_admm_charge_constrained_projection
     369            2 :             dft_control%admm_control%scaling_model = do_admm_exch_scaling_none
     370              :          CASE DEFAULT
     371              :             CALL cp_abort(__LOCATION__, &
     372          524 :                           "ADMM_TYPE keyword in &AUXILIARY_DENSITY_MATRIX_METHOD section has a wrong value.")
     373              :          END SELECT
     374              : 
     375              :          CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%EPS_FILTER", &
     376          524 :                                    r_val=dft_control%admm_control%eps_filter)
     377              : 
     378          524 :          CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%EXCH_CORRECTION_FUNC", i_val=method_id)
     379          524 :          dft_control%admm_control%aux_exch_func = method_id
     380              : 
     381              :          ! parameters for X functional
     382          524 :          dft_control%admm_control%aux_exch_func_param = .FALSE.
     383              :          CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%OPTX_A1", explicit=explicit, &
     384          524 :                                    r_val=dft_control%admm_control%aux_x_param(1))
     385          524 :          IF (explicit) dft_control%admm_control%aux_exch_func_param = .TRUE.
     386              :          CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%OPTX_A2", explicit=explicit, &
     387          524 :                                    r_val=dft_control%admm_control%aux_x_param(2))
     388          524 :          IF (explicit) dft_control%admm_control%aux_exch_func_param = .TRUE.
     389              :          CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%OPTX_GAMMA", explicit=explicit, &
     390          524 :                                    r_val=dft_control%admm_control%aux_x_param(3))
     391          524 :          IF (explicit) dft_control%admm_control%aux_exch_func_param = .TRUE.
     392              : 
     393          524 :          CALL read_admm_block_list(dft_control%admm_control, dft_section)
     394              : 
     395              :          ! check for double assignments
     396            2 :          SELECT CASE (admmtype)
     397              :          CASE (admm2_type)
     398            2 :             IF (exopt2) CALL cp_warn(__LOCATION__, &
     399            0 :                                      "Value of ADMM_PURIFICATION_METHOD keyword will be overwritten with ADMM_TYPE selections.")
     400            2 :             IF (exopt3) CALL cp_warn(__LOCATION__, &
     401            0 :                                      "Value of EXCH_SCALING_MODEL keyword will be overwritten with ADMM_TYPE selections.")
     402              :          CASE (admm1_type, admms_type, admmp_type, admmq_type)
     403           16 :             IF (exopt1) CALL cp_warn(__LOCATION__, &
     404            2 :                                      "Value of METHOD keyword will be overwritten with ADMM_TYPE selections.")
     405           16 :             IF (exopt2) CALL cp_warn(__LOCATION__, &
     406            2 :                                      "Value of METHOD keyword will be overwritten with ADMM_TYPE selections.")
     407           16 :             IF (exopt3) CALL cp_warn(__LOCATION__, &
     408          524 :                                      "Value of EXCH_SCALING_MODEL keyword will be overwritten with ADMM_TYPE selections.")
     409              :          END SELECT
     410              : 
     411              :          !    In the case of charge-constrained projection (e.g. according to Merlot),
     412              :          !    there is no purification needed and hence, do_admm_purify_none has to be set.
     413              : 
     414              :          IF ((dft_control%admm_control%method == do_admm_blocking_purify_full .OR. &
     415              :               dft_control%admm_control%method == do_admm_blocked_projection) &
     416          524 :              .AND. dft_control%admm_control%scaling_model == do_admm_exch_scaling_merlot) THEN
     417            0 :             CPABORT("ADMM: Blocking and Merlot scaling are mutually exclusive.")
     418              :          END IF
     419              : 
     420          524 :          IF (dft_control%admm_control%method == do_admm_charge_constrained_projection .AND. &
     421              :              dft_control%admm_control%purification_method /= do_admm_purify_none) THEN
     422              :             CALL cp_abort(__LOCATION__, &
     423              :                           "ADMM: In the case of METHOD=CHARGE_CONSTRAINED_PROJECTION, "// &
     424            0 :                           "ADMM_PURIFICATION_METHOD=NONE has to be set.")
     425              :          END IF
     426              : 
     427          524 :          IF (dft_control%admm_control%purification_method == do_admm_purify_mo_diag .OR. &
     428              :              dft_control%admm_control%purification_method == do_admm_purify_mo_no_diag) THEN
     429           62 :             IF (dft_control%admm_control%method /= do_admm_basis_projection) THEN
     430            0 :                CPABORT("ADMM: Chosen purification requires BASIS_PROJECTION")
     431              :             END IF
     432              : 
     433           62 :             IF (.NOT. do_ot) CPABORT("ADMM: MO-based purification requires OT.")
     434              :          END IF
     435              : 
     436          524 :          IF (dft_control%admm_control%purification_method == do_admm_purify_none_dm .OR. &
     437              :              dft_control%admm_control%purification_method == do_admm_purify_mcweeny) THEN
     438           14 :             dft_control%do_admm_dm = .TRUE.
     439              :          ELSE
     440          510 :             dft_control%do_admm_mo = .TRUE.
     441              :          END IF
     442              :       END IF
     443              : 
     444              :       ! Set restricted to true, if both OT and ROKS are requested
     445              :       !MK in principle dft_control%restricted could be dropped completely like the
     446              :       !MK input key by using only dft_control%roks now
     447         9032 :       CALL section_vals_val_get(scf_section, "OT%_SECTION_PARAMETERS_", l_val=l_param)
     448         9032 :       dft_control%restricted = (dft_control%roks .AND. l_param)
     449              : 
     450         9032 :       CALL section_vals_val_get(dft_section, "CHARGE", i_val=dft_control%charge)
     451         9032 :       CALL section_vals_val_get(dft_section, "MULTIPLICITY", i_val=dft_control%multiplicity)
     452         9032 :       CALL section_vals_val_get(dft_section, "RELAX_MULTIPLICITY", r_val=dft_control%relax_multiplicity)
     453         9032 :       IF (dft_control%relax_multiplicity > 0.0_dp) THEN
     454           10 :          IF (.NOT. dft_control%uks) THEN
     455              :             CALL cp_abort(__LOCATION__, "The option RELAX_MULTIPLICITY is only valid for "// &
     456            0 :                           "unrestricted Kohn-Sham (UKS) calculations")
     457              :          END IF
     458              :       END IF
     459              : 
     460              :       !Read the HAIR PROBES input section if present
     461         9032 :       hairy_probes_section => section_vals_get_subs_vals(dft_section, "HAIRY_PROBES")
     462         9032 :       CALL section_vals_get(hairy_probes_section, n_repetition=nrep, explicit=is_present)
     463              : 
     464         9032 :       IF (is_present) THEN
     465            4 :          dft_control%hairy_probes = .TRUE.
     466           20 :          ALLOCATE (dft_control%probe(nrep))
     467            4 :          CALL read_hairy_probes_sections(dft_control, hairy_probes_section)
     468              :       END IF
     469              : 
     470              :       ! check for the presence of the low spin roks section
     471         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "LOW_SPIN_ROKS")
     472         9032 :       CALL section_vals_get(tmp_section, explicit=dft_control%low_spin_roks)
     473              : 
     474         9032 :       dft_control%sic_method_id = sic_none
     475         9032 :       dft_control%sic_scaling_a = 1.0_dp
     476         9032 :       dft_control%sic_scaling_b = 1.0_dp
     477              : 
     478              :       ! DFT+U
     479         9032 :       dft_control%dft_plus_u = .FALSE.
     480         9032 :       CALL section_vals_val_get(dft_section, "PLUS_U_METHOD", i_val=method_id)
     481         9032 :       dft_control%plus_u_method_id = method_id
     482              : 
     483              :       ! Minimum tracking linear response U and J
     484         9032 :       CALL section_vals_val_get(dft_section, "EPS_U_J_LOOP", r_val=dft_control%eps_u_j_loop)
     485         9032 :       CALL section_vals_val_get(dft_section, "MAX_MTLR_LOOP", i_val=dft_control%max_mtlr_iter)
     486              : 
     487              :       ! Smearing in use
     488         9032 :       dft_control%smear = .FALSE.
     489              : 
     490              :       ! Surface dipole correction
     491         9032 :       dft_control%correct_surf_dip = .FALSE.
     492         9032 :       CALL section_vals_val_get(dft_section, "SURFACE_DIPOLE_CORRECTION", l_val=dft_control%correct_surf_dip)
     493         9032 :       CALL section_vals_val_get(dft_section, "SURF_DIP_DIR", i_val=dft_control%dir_surf_dip)
     494         9032 :       dft_control%pos_dir_surf_dip = -1.0_dp
     495         9032 :       CALL section_vals_val_get(dft_section, "SURF_DIP_POS", r_val=dft_control%pos_dir_surf_dip)
     496              :       ! another logical variable, surf_dip_correct_switch, is introduced for
     497              :       ! implementation of "SURF_DIP_SWITCH" [SGh]
     498         9032 :       dft_control%switch_surf_dip = .FALSE.
     499         9032 :       dft_control%surf_dip_correct_switch = dft_control%correct_surf_dip
     500         9032 :       CALL section_vals_val_get(dft_section, "SURF_DIP_SWITCH", l_val=dft_control%switch_surf_dip)
     501         9032 :       dft_control%correct_el_density_dip = .FALSE.
     502         9032 :       CALL section_vals_val_get(dft_section, "CORE_CORR_DIP", l_val=dft_control%correct_el_density_dip)
     503         9032 :       IF (dft_control%correct_el_density_dip) THEN
     504            4 :          IF (dft_control%correct_surf_dip) THEN
     505              :             ! Do nothing, move on
     506              :          ELSE
     507            0 :             dft_control%correct_el_density_dip = .FALSE.
     508            0 :             CPWARN("CORE_CORR_DIP keyword is activated only if SURFACE_DIPOLE_CORRECTION is TRUE")
     509              :          END IF
     510              :       END IF
     511              : 
     512              :       CALL section_vals_val_get(dft_section, "BASIS_SET_FILE_NAME", &
     513         9032 :                                 c_val=basis_set_file_name)
     514              :       CALL section_vals_val_get(dft_section, "POTENTIAL_FILE_NAME", &
     515         9032 :                                 c_val=potential_file_name)
     516              : 
     517              :       ! Read the input section
     518         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "sic")
     519              :       CALL section_vals_val_get(tmp_section, "SIC_METHOD", &
     520         9032 :                                 i_val=dft_control%sic_method_id)
     521              :       CALL section_vals_val_get(tmp_section, "ORBITAL_SET", &
     522         9032 :                                 i_val=dft_control%sic_list_id)
     523              :       CALL section_vals_val_get(tmp_section, "SIC_SCALING_A", &
     524         9032 :                                 r_val=dft_control%sic_scaling_a)
     525              :       CALL section_vals_val_get(tmp_section, "SIC_SCALING_B", &
     526         9032 :                                 r_val=dft_control%sic_scaling_b)
     527              : 
     528         9032 :       do_rtp = .FALSE.
     529         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "REAL_TIME_PROPAGATION")
     530         9032 :       CALL section_vals_get(tmp_section, explicit=is_present)
     531         9032 :       IF (is_present) THEN
     532          324 :          CALL read_rtp_section(dft_control, tmp_section)
     533          324 :          do_rtp = .TRUE.
     534              :       END IF
     535              : 
     536              :       ! Read the input section
     537         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "XAS")
     538         9032 :       CALL section_vals_get(tmp_section, explicit=dft_control%do_xas_calculation)
     539         9032 :       IF (dft_control%do_xas_calculation) THEN
     540              :          ! Override with section parameter
     541              :          CALL section_vals_val_get(tmp_section, "_SECTION_PARAMETERS_", &
     542           42 :                                    l_val=dft_control%do_xas_calculation)
     543              :       END IF
     544              : 
     545         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "XAS_TDP")
     546         9032 :       CALL section_vals_get(tmp_section, explicit=dft_control%do_xas_tdp_calculation)
     547         9032 :       IF (dft_control%do_xas_tdp_calculation) THEN
     548              :          ! Override with section parameter
     549              :          CALL section_vals_val_get(tmp_section, "_SECTION_PARAMETERS_", &
     550           52 :                                    l_val=dft_control%do_xas_tdp_calculation)
     551              :       END IF
     552              : 
     553              :       ! Read the finite field input section
     554         9032 :       dft_control%apply_efield = .FALSE.
     555         9032 :       dft_control%apply_efield_field = .FALSE. !this is for RTP
     556         9032 :       dft_control%apply_vector_potential = .FALSE. !this is for RTP
     557         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "EFIELD")
     558         9032 :       CALL section_vals_get(tmp_section, n_repetition=nrep, explicit=is_present)
     559         9032 :       IF (is_present) THEN
     560         1368 :          ALLOCATE (dft_control%efield_fields(nrep))
     561          342 :          CALL read_efield_sections(dft_control, tmp_section, cell)
     562          342 :          IF (do_rtp) THEN
     563           30 :             IF (.NOT. dft_control%rtp_control%velocity_gauge) THEN
     564           20 :                dft_control%apply_efield_field = .TRUE.
     565              :             ELSE
     566           10 :                dft_control%apply_vector_potential = .TRUE.
     567              :                ! Use this input value of vector potential to (re)start RTP
     568           40 :                dft_control%rtp_control%vec_pot = dft_control%efield_fields(1)%efield%vec_pot_initial
     569              :             END IF
     570              :          ELSE
     571          312 :             dft_control%apply_efield = .TRUE.
     572          312 :             CPASSERT(nrep == 1)
     573              :          END IF
     574              :       END IF
     575              : 
     576              :       ! Now, can try to guess polarisation in rtp
     577         9032 :       IF (do_rtp) THEN
     578              :          ! tmp_section => section_vals_get_subs_vals(dft_section, "REAL_TIME_PROPAGATION%PRINT%POLARIZABILITY")
     579              :          ! CALL section_vals_get(tmp_section, explicit=is_present)
     580              :          local_moment_possible = (dft_control%rtp_control%rtp_method == rtp_method_bse .OR. &
     581              :                                   dft_control%rtp_control%rtp_method == rtp_method_bse_linearized) .OR. &
     582          324 :                                  ((.NOT. dft_control%rtp_control%periodic) .AND. dft_control%rtp_control%linear_scaling)
     583           90 :          IF (local_moment_possible .AND. (.NOT. ASSOCIATED(dft_control%rtp_control%print_pol_elements))) THEN
     584           90 :             tmp_section => section_vals_get_subs_vals(dft_section, "REAL_TIME_PROPAGATION")
     585              :             CALL guess_pol_elements(dft_control, &
     586           90 :                                     dft_control%rtp_control%print_pol_elements)
     587              :          END IF
     588              :       END IF
     589              : 
     590              :       ! Read the finite field input section for periodic fields
     591         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "PERIODIC_EFIELD")
     592         9032 :       CALL section_vals_get(tmp_section, explicit=dft_control%apply_period_efield)
     593         9032 :       IF (dft_control%apply_period_efield) THEN
     594          532 :          ALLOCATE (dft_control%period_efield)
     595           76 :          CALL section_vals_val_get(tmp_section, "POLARISATION", r_vals=pol)
     596          532 :          dft_control%period_efield%polarisation(1:3) = pol(1:3)
     597           76 :          IF (PRESENT(cell)) THEN
     598           76 :             IF (ASSOCIATED(cell)) THEN
     599           76 :                CALL cell_transform_input_cartesian(cell, dft_control%period_efield%polarisation(1:3))
     600              :             END IF
     601              :          END IF
     602           76 :          CALL section_vals_val_get(tmp_section, "D_FILTER", r_vals=pol)
     603          532 :          dft_control%period_efield%d_filter(1:3) = pol(1:3)
     604           76 :          IF (PRESENT(cell)) THEN
     605           76 :             IF (ASSOCIATED(cell)) CALL cell_transform_input_cartesian(cell, dft_control%period_efield%d_filter(1:3))
     606              :          END IF
     607              :          CALL section_vals_val_get(tmp_section, "INTENSITY", &
     608           76 :                                    r_val=dft_control%period_efield%strength)
     609           76 :          dft_control%period_efield%displacement_field = .FALSE.
     610              :          CALL section_vals_val_get(tmp_section, "DISPLACEMENT_FIELD", &
     611           76 :                                    l_val=dft_control%period_efield%displacement_field)
     612              : 
     613           76 :          CALL section_vals_val_get(tmp_section, "INTENSITY_LIST", r_vals=pol)
     614              : 
     615           76 :          CALL section_vals_val_get(tmp_section, "INTENSITIES_FILE_NAME", c_val=intensities_file_name)
     616              : 
     617           76 :          IF (SIZE(pol) > 1 .OR. pol(1) /= 0.0_dp) THEN
     618              :             ! if INTENSITY_LIST is present, INTENSITY and INTENSITIES_FILE_NAME must not be present
     619            2 :             IF (dft_control%period_efield%strength /= 0.0_dp .OR. intensities_file_name /= "") THEN
     620              :                CALL cp_abort(__LOCATION__, "[PERIODIC FIELD] Only one of INTENSITY, INTENSITY_LIST "// &
     621            0 :                              "or INTENSITIES_FILE_NAME can be specified.")
     622              :             END IF
     623              : 
     624            6 :             ALLOCATE (dft_control%period_efield%strength_list(SIZE(pol)))
     625           50 :             dft_control%period_efield%strength_list(1:SIZE(pol)) = pol(1:SIZE(pol))
     626              :          END IF
     627              : 
     628           76 :          IF (intensities_file_name /= "") THEN
     629              :             ! if INTENSITIES_FILE_NAME is present, INTENSITY must not be present
     630            2 :             IF (dft_control%period_efield%strength /= 0.0_dp) THEN
     631              :                CALL cp_abort(__LOCATION__, "[PERIODIC FIELD] Only one of INTENSITY, INTENSITY_LIST "// &
     632            0 :                              "or INTENSITIES_FILE_NAME can be specified.")
     633              :             END IF
     634              : 
     635            2 :             CALL parser_create(parser, intensities_file_name)
     636              : 
     637            2 :             nrep = 0
     638           24 :             DO WHILE (.TRUE.)
     639           26 :                CALL parser_read_line(parser, 1, at_end)
     640           26 :                IF (at_end) EXIT
     641           24 :                nrep = nrep + 1
     642              :             END DO
     643              : 
     644            2 :             IF (nrep == 0) THEN
     645            0 :                CPABORT("[PERIODIC FIELD] No intensities found in INTENSITIES_FILE_NAME")
     646              :             END IF
     647              : 
     648            6 :             ALLOCATE (dft_control%period_efield%strength_list(nrep))
     649              : 
     650            2 :             CALL parser_reset(parser)
     651           26 :             DO irep = 1, nrep
     652           24 :                CALL parser_read_line(parser, 1)
     653           26 :                READ (parser%input_line, *) dft_control%period_efield%strength_list(irep)
     654              :             END DO
     655              : 
     656            4 :             CALL parser_release(parser)
     657              :          END IF
     658              : 
     659              :          CALL section_vals_val_get(tmp_section, "START_FRAME", &
     660           76 :                                    i_val=dft_control%period_efield%start_frame)
     661              :          CALL section_vals_val_get(tmp_section, "END_FRAME", &
     662           76 :                                    i_val=dft_control%period_efield%end_frame)
     663              : 
     664           76 :          IF (dft_control%period_efield%end_frame /= -1) THEN
     665              :             ! check if valid bounds are given
     666              :             ! if an end frame is given, the number of active frames must be a
     667              :             ! multiple of the number of intensities
     668            4 :             IF (dft_control%period_efield%start_frame > dft_control%period_efield%end_frame) THEN
     669            0 :                CPABORT("[PERIODIC FIELD] START_FRAME > END_FRAME")
     670            4 :             ELSE IF (dft_control%period_efield%start_frame < 1) THEN
     671            0 :                CPABORT("[PERIODIC FIELD] START_FRAME < 1")
     672            4 :             ELSE IF (MOD(dft_control%period_efield%end_frame - &
     673              :                          dft_control%period_efield%start_frame + 1, SIZE(pol)) /= 0) THEN
     674              :                CALL cp_abort(__LOCATION__, &
     675            0 :                              "[PERIODIC FIELD] Number of active frames must be a multiple of the number of intensities")
     676              :             END IF
     677              :          END IF
     678              : 
     679              :          ! periodic fields don't work with RTP
     680           76 :          IF (do_rtp) THEN
     681              :             CALL cp_abort(__LOCATION__, &
     682              :                           "Periodic efield cannot be used with RTP. When restarting a "// &
     683              :                           "run with periodic efield, set RESTART_RTP under &EXT_RESTART "// &
     684            0 :                           "section to .FALSE. explicitly if RESTART_DEFAULT is .TRUE.")
     685              :          END IF
     686           76 :          IF (dft_control%period_efield%displacement_field) THEN
     687           16 :             CALL cite_reference(Stengel2009)
     688              :          ELSE
     689           60 :             CALL cite_reference(Souza2002)
     690           60 :             CALL cite_reference(Umari2002)
     691              :          END IF
     692              :       END IF
     693              : 
     694              :       ! Read the external potential input section
     695         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "EXTERNAL_POTENTIAL")
     696         9032 :       CALL section_vals_get(tmp_section, explicit=dft_control%apply_external_potential)
     697         9032 :       IF (dft_control%apply_external_potential) THEN
     698           16 :          CALL expot_control_create(dft_control%expot_control)
     699              :          CALL section_vals_val_get(tmp_section, "READ_FROM_CUBE", &
     700           16 :                                    l_val=dft_control%expot_control%read_from_cube)
     701              :          CALL section_vals_val_get(tmp_section, "STATIC", &
     702           16 :                                    l_val=dft_control%expot_control%static)
     703              :          CALL section_vals_val_get(tmp_section, "SCALING_FACTOR", &
     704           16 :                                    r_val=dft_control%expot_control%scaling_factor)
     705              :          ! External potential using Maxwell equation
     706           16 :          maxwell_section => section_vals_get_subs_vals(tmp_section, "MAXWELL")
     707           16 :          CALL section_vals_get(maxwell_section, explicit=is_present)
     708           16 :          IF (is_present) THEN
     709            0 :             dft_control%expot_control%maxwell_solver = .TRUE.
     710            0 :             CALL maxwell_control_create(dft_control%maxwell_control)
     711              :             ! read the input values from Maxwell section
     712              :             CALL section_vals_val_get(maxwell_section, "TEST_REAL", &
     713            0 :                                       r_val=dft_control%maxwell_control%real_test)
     714              :             CALL section_vals_val_get(maxwell_section, "TEST_INTEGER", &
     715            0 :                                       i_val=dft_control%maxwell_control%int_test)
     716              :             CALL section_vals_val_get(maxwell_section, "TEST_LOGICAL", &
     717            0 :                                       l_val=dft_control%maxwell_control%log_test)
     718              :          ELSE
     719           16 :             dft_control%expot_control%maxwell_solver = .FALSE.
     720              :          END IF
     721              :       END IF
     722              : 
     723              :       ! Read the SCCS input section if present
     724         9032 :       sccs_section => section_vals_get_subs_vals(dft_section, "SCCS")
     725         9032 :       CALL section_vals_get(sccs_section, explicit=is_present)
     726         9032 :       IF (is_present) THEN
     727              :          ! Check section parameter if SCCS is activated
     728              :          CALL section_vals_val_get(sccs_section, "_SECTION_PARAMETERS_", &
     729           12 :                                    l_val=dft_control%do_sccs)
     730           12 :          IF (dft_control%do_sccs) THEN
     731           12 :             ALLOCATE (dft_control%sccs_control)
     732              :             CALL section_vals_val_get(sccs_section, "RELATIVE_PERMITTIVITY", &
     733           12 :                                       r_val=dft_control%sccs_control%epsilon_solvent)
     734              :             CALL section_vals_val_get(sccs_section, "ALPHA", &
     735           12 :                                       r_val=dft_control%sccs_control%alpha_solvent)
     736              :             CALL section_vals_val_get(sccs_section, "BETA", &
     737           12 :                                       r_val=dft_control%sccs_control%beta_solvent)
     738              :             CALL section_vals_val_get(sccs_section, "DELTA_RHO", &
     739           12 :                                       r_val=dft_control%sccs_control%delta_rho)
     740              :             CALL section_vals_val_get(sccs_section, "DERIVATIVE_METHOD", &
     741           12 :                                       i_val=dft_control%sccs_control%derivative_method)
     742              :             CALL section_vals_val_get(sccs_section, "METHOD", &
     743           12 :                                       i_val=dft_control%sccs_control%method_id)
     744              :             CALL section_vals_val_get(sccs_section, "EPS_SCCS", &
     745           12 :                                       r_val=dft_control%sccs_control%eps_sccs)
     746              :             CALL section_vals_val_get(sccs_section, "EPS_SCF", &
     747           12 :                                       r_val=dft_control%sccs_control%eps_scf)
     748              :             CALL section_vals_val_get(sccs_section, "GAMMA", &
     749           12 :                                       r_val=dft_control%sccs_control%gamma_solvent)
     750              :             CALL section_vals_val_get(sccs_section, "MAX_ITER", &
     751           12 :                                       i_val=dft_control%sccs_control%max_iter)
     752              :             CALL section_vals_val_get(sccs_section, "MIXING", &
     753           12 :                                       r_val=dft_control%sccs_control%mixing)
     754           22 :             SELECT CASE (dft_control%sccs_control%method_id)
     755              :             CASE (sccs_andreussi)
     756           10 :                tmp_section => section_vals_get_subs_vals(sccs_section, "ANDREUSSI")
     757              :                CALL section_vals_val_get(tmp_section, "RHO_MAX", &
     758           10 :                                          r_val=dft_control%sccs_control%rho_max)
     759              :                CALL section_vals_val_get(tmp_section, "RHO_MIN", &
     760           10 :                                          r_val=dft_control%sccs_control%rho_min)
     761           10 :                IF (dft_control%sccs_control%rho_max < dft_control%sccs_control%rho_min) THEN
     762              :                   CALL cp_abort(__LOCATION__, &
     763              :                                 "The SCCS parameter RHO_MAX is smaller than RHO_MIN. "// &
     764            0 :                                 "Please, check your input!")
     765              :                END IF
     766           10 :                CALL cite_reference(Andreussi2012)
     767              :             CASE (sccs_fattebert_gygi)
     768            2 :                tmp_section => section_vals_get_subs_vals(sccs_section, "FATTEBERT-GYGI")
     769              :                CALL section_vals_val_get(tmp_section, "BETA", &
     770            2 :                                          r_val=dft_control%sccs_control%beta)
     771            2 :                IF (dft_control%sccs_control%beta < 0.5_dp) THEN
     772              :                   CALL cp_abort(__LOCATION__, &
     773              :                                 "A value smaller than 0.5 for the SCCS parameter beta "// &
     774            0 :                                 "causes numerical problems. Please, check your input!")
     775              :                END IF
     776              :                CALL section_vals_val_get(tmp_section, "RHO_ZERO", &
     777            2 :                                          r_val=dft_control%sccs_control%rho_zero)
     778            2 :                CALL cite_reference(Fattebert2002)
     779              :             CASE (sccs_saa_andreussi)
     780            0 :                tmp_section => section_vals_get_subs_vals(sccs_section, "SAA_ANDREUSSI")
     781            0 :                CALL section_vals_get(tmp_section, explicit=is_present)
     782            0 :                IF (.NOT. is_present) THEN
     783              :                   CALL cp_abort(__LOCATION__, &
     784              :                                 "SCCS method SAA_ANDREUSSI requires the "// &
     785            0 :                                 "SAA_ANDREUSSI section.")
     786              :                END IF
     787            0 :                IF (.NOT. ALL(cell%perd == 1)) THEN
     788              :                   CALL cp_abort(__LOCATION__, &
     789              :                                 "SCCS method SAA_ANDREUSSI is only implemented for "// &
     790            0 :                                 "3D periodic calculations.")
     791              :                END IF
     792              :                CALL section_vals_val_get(tmp_section, "RHO_MAX", &
     793            0 :                                          r_val=dft_control%sccs_control%rho_max)
     794              :                CALL section_vals_val_get(tmp_section, "RHO_MIN", &
     795            0 :                                          r_val=dft_control%sccs_control%rho_min)
     796            0 :                IF (dft_control%sccs_control%rho_max < dft_control%sccs_control%rho_min) THEN
     797              :                   CALL cp_abort(__LOCATION__, &
     798              :                                 "The SCCS parameter RHO_MAX is smaller than RHO_MIN. "// &
     799            0 :                                 "Please, check your input!")
     800              :                END IF
     801              :                CALL section_vals_val_get(tmp_section, "F0", &
     802            0 :                                          r_val=dft_control%sccs_control%f0)
     803              :                CALL section_vals_val_get(tmp_section, "DELTA_ETA", &
     804            0 :                                          r_val=dft_control%sccs_control%delta_eta)
     805              :                CALL section_vals_val_get(tmp_section, "DELTA_ZETA", &
     806            0 :                                          r_val=dft_control%sccs_control%delta_zeta)
     807              :                CALL section_vals_val_get(tmp_section, "R_SOLV", &
     808            0 :                                          r_val=dft_control%sccs_control%r_solv)
     809              :                CALL section_vals_val_get(tmp_section, "ALPHA_ZETA", &
     810            0 :                                          r_val=dft_control%sccs_control%alpha_zeta)
     811            0 :                CALL cite_reference(Andreussi2019)
     812            0 :                CALL cite_reference(Chai2025a)
     813              :             CASE DEFAULT
     814           12 :                CPABORT("Invalid SCCS model specified. Please, check your input!")
     815              :             END SELECT
     816           12 :             CALL cite_reference(Yin2017)
     817              :          END IF
     818              :       END IF
     819              : 
     820              :       ! Read the planar counter charge section
     821         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "PLANAR_COUNTER_CHARGE")
     822         9032 :       CALL section_vals_get(tmp_section, explicit=is_present)
     823         9032 :       IF (is_present) THEN
     824              :          ! Check section parameter if planar counter charge is activated
     825              :          CALL section_vals_val_get(tmp_section, "_SECTION_PARAMETERS_", &
     826            6 :                                    l_val=dft_control%do_pcc)
     827            6 :          IF (dft_control%do_pcc) THEN
     828            6 :             ALLOCATE (dft_control%pcc_control)
     829              :             CALL section_vals_val_get(tmp_section, "DIST_EDGE", &
     830            6 :                                       r_val=dft_control%pcc_control%dist_edge)
     831              :             CALL section_vals_val_get(tmp_section, "GAU_C", &
     832            6 :                                       r_val=dft_control%pcc_control%gau_c)
     833              :             CALL section_vals_val_get(tmp_section, "PARALLEL_PLANE", &
     834            6 :                                       i_val=dft_control%pcc_control%surf_normal)
     835              :          END IF
     836              :       END IF
     837              : 
     838              :       ! Read the planar averaged Hartree potential section
     839         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "PLANAR_AVERAGED_V_HARTREE")
     840         9032 :       CALL section_vals_get(tmp_section, explicit=is_present)
     841         9032 :       IF (is_present) THEN
     842              :          ! Check section parameter if the planar-averaged potential is activated
     843              :          CALL section_vals_val_get(tmp_section, "_SECTION_PARAMETERS_", &
     844            2 :                                    l_val=dft_control%do_paep)
     845            2 :          IF (dft_control%do_paep) THEN
     846            2 :             ALLOCATE (dft_control%paep_control)
     847              :             CALL section_vals_val_get(tmp_section, "PARALLEL_PLANE", &
     848            2 :                                       i_val=dft_control%paep_control%surf_normal)
     849              :          END IF
     850              :       END IF
     851              : 
     852              :       ! ZMP added input sections
     853              :       ! Read the external density input section
     854         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "EXTERNAL_DENSITY")
     855         9032 :       CALL section_vals_get(tmp_section, explicit=dft_control%apply_external_density)
     856              : 
     857              :       ! Read the external vxc input section
     858         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "EXTERNAL_VXC")
     859         9032 :       CALL section_vals_get(tmp_section, explicit=dft_control%apply_external_vxc)
     860              : 
     861         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "LOCALIZE")
     862         9032 :       CALL section_vals_val_get(tmp_section, "EACH", i_val=dft_control%localize_each)
     863              : 
     864              :       ! SMEAGOL interface
     865         9032 :       tmp_section => section_vals_get_subs_vals(dft_section, "SMEAGOL")
     866         9032 :       CALL read_smeagol_control(dft_control%smeagol_control, tmp_section)
     867              : 
     868        27096 :    END SUBROUTINE read_dft_control
     869              : 
     870              : ! **************************************************************************************************
     871              : !> \brief Reads the input and stores in the rixs_control_type
     872              : !> \param rixs_control ...
     873              : !> \param rixs_section ...
     874              : !> \param qs_control ...
     875              : ! **************************************************************************************************
     876           32 :    SUBROUTINE read_rixs_control(rixs_control, rixs_section, qs_control)
     877              :       TYPE(rixs_control_type), POINTER                   :: rixs_control
     878              :       TYPE(section_vals_type), POINTER                   :: rixs_section
     879              :       TYPE(qs_control_type), POINTER                     :: qs_control
     880              : 
     881              :       TYPE(section_vals_type), POINTER                   :: td_section, xas_section
     882              : 
     883           32 :       CALL section_vals_val_get(rixs_section, "_SECTION_PARAMETERS_", l_val=rixs_control%enabled)
     884              : 
     885           32 :       CALL section_vals_val_get(rixs_section, "CORE_STATES", i_val=rixs_control%core_states)
     886           32 :       CALL section_vals_val_get(rixs_section, "VALENCE_STATES", i_val=rixs_control%valence_states)
     887              : 
     888           32 :       td_section => section_vals_get_subs_vals(rixs_section, "TDDFPT")
     889           32 :       CALL read_tddfpt2_control(rixs_control%tddfpt2_control, td_section, qs_control)
     890              : 
     891           32 :       xas_section => section_vals_get_subs_vals(rixs_section, "XAS_TDP")
     892           32 :       CALL read_xas_tdp_control(rixs_control%xas_tdp_control, xas_section)
     893              : 
     894           32 :    END SUBROUTINE read_rixs_control
     895              : 
     896              : ! **************************************************************************************************
     897              : !> \brief ...
     898              : !> \param qs_control ...
     899              : !> \param dft_section ...
     900              : ! **************************************************************************************************
     901         9032 :    SUBROUTINE read_mgrid_section(qs_control, dft_section)
     902              : 
     903              :       TYPE(qs_control_type), INTENT(INOUT)               :: qs_control
     904              :       TYPE(section_vals_type), POINTER                   :: dft_section
     905              : 
     906              :       CHARACTER(len=*), PARAMETER :: routineN = 'read_mgrid_section'
     907              : 
     908              :       INTEGER                                            :: handle, igrid_level, interp_kind, &
     909              :                                                             ngrid_level
     910              :       LOGICAL                                            :: explicit, multigrid_set
     911              :       REAL(dp)                                           :: cutoff
     912         9032 :       REAL(dp), DIMENSION(:), POINTER                    :: cutofflist
     913              :       TYPE(section_vals_type), POINTER                   :: interp_section, mgrid_section
     914              : 
     915         9032 :       CALL timeset(routineN, handle)
     916              : 
     917         9032 :       NULLIFY (interp_section, mgrid_section, cutofflist)
     918         9032 :       mgrid_section => section_vals_get_subs_vals(dft_section, "MGRID")
     919         9032 :       interp_section => section_vals_get_subs_vals(mgrid_section, "INTERPOLATOR")
     920              : 
     921         9032 :       CALL section_vals_val_get(mgrid_section, "NGRIDS", i_val=ngrid_level)
     922         9032 :       CALL section_vals_val_get(mgrid_section, "MULTIGRID_SET", l_val=multigrid_set)
     923         9032 :       CALL section_vals_val_get(mgrid_section, "CUTOFF", r_val=cutoff)
     924         9032 :       CALL section_vals_val_get(mgrid_section, "PROGRESSION_FACTOR", r_val=qs_control%progression_factor)
     925         9032 :       CALL section_vals_val_get(mgrid_section, "COMMENSURATE", l_val=qs_control%commensurate_mgrids)
     926         9032 :       CALL section_vals_val_get(interp_section, "KIND", i_val=interp_kind)
     927         9032 :       IF (interp_kind /= pw_interp) qs_control%commensurate_mgrids = .TRUE.
     928         9032 :       CALL section_vals_val_get(mgrid_section, "REALSPACE", l_val=qs_control%realspace_mgrids)
     929         9032 :       CALL section_vals_val_get(mgrid_section, "REL_CUTOFF", r_val=qs_control%relative_cutoff)
     930              :       CALL section_vals_val_get(mgrid_section, "SKIP_LOAD_BALANCE_DISTRIBUTED", &
     931         9032 :                                 l_val=qs_control%skip_load_balance_distributed)
     932              : 
     933              :       ! For SE and DFTB possibly override with new defaults
     934         9032 :       IF (qs_control%semi_empirical .OR. qs_control%dftb .OR. qs_control%xtb) THEN
     935         2526 :          ngrid_level = 1
     936         2526 :          multigrid_set = .FALSE.
     937              :          ! Override default cutoff value unless user specified an explicit argument..
     938         2526 :          CALL section_vals_val_get(mgrid_section, "CUTOFF", explicit=explicit, r_val=cutoff)
     939         2526 :          IF (.NOT. explicit) cutoff = 1.0_dp
     940              :       END IF
     941              : 
     942        27096 :       ALLOCATE (qs_control%e_cutoff(ngrid_level))
     943         9032 :       qs_control%cutoff = cutoff
     944              : 
     945         9032 :       IF (multigrid_set) THEN
     946              :          ! Read the values from input
     947            4 :          IF (qs_control%commensurate_mgrids) THEN
     948            0 :             CPABORT("Do not specify cutoffs for the commensurate grids (NYI)")
     949              :          END IF
     950              : 
     951            4 :          CALL section_vals_val_get(mgrid_section, "MULTIGRID_CUTOFF", r_vals=cutofflist)
     952            4 :          IF (ASSOCIATED(cutofflist)) THEN
     953            4 :             IF (SIZE(cutofflist, 1) /= ngrid_level) THEN
     954            0 :                CPABORT("Number of multi-grids requested and number of cutoff values do not match")
     955              :             END IF
     956           20 :             DO igrid_level = 1, ngrid_level
     957           20 :                qs_control%e_cutoff(igrid_level) = cutofflist(igrid_level)
     958              :             END DO
     959              :          END IF
     960              :          ! set cutoff to smallest value in multgrid available with >= cutoff
     961           20 :          DO igrid_level = ngrid_level, 1, -1
     962           16 :             IF (qs_control%cutoff <= qs_control%e_cutoff(igrid_level)) THEN
     963            0 :                qs_control%cutoff = qs_control%e_cutoff(igrid_level)
     964            0 :                EXIT
     965              :             END IF
     966              :             ! set largest grid value to cutoff
     967           20 :             IF (igrid_level == 1) THEN
     968            4 :                qs_control%cutoff = qs_control%e_cutoff(1)
     969              :             END IF
     970              :          END DO
     971              :       ELSE
     972         9028 :          IF (qs_control%commensurate_mgrids) qs_control%progression_factor = 4.0_dp
     973         9028 :          qs_control%e_cutoff(1) = qs_control%cutoff
     974        28250 :          DO igrid_level = 2, ngrid_level
     975              :             qs_control%e_cutoff(igrid_level) = qs_control%e_cutoff(igrid_level - 1)/ &
     976        28250 :                                                qs_control%progression_factor
     977              :          END DO
     978              :       END IF
     979              :       ! check that multigrids are ordered
     980        28266 :       DO igrid_level = 2, ngrid_level
     981        28266 :          IF (qs_control%e_cutoff(igrid_level) > qs_control%e_cutoff(igrid_level - 1)) THEN
     982            0 :             CPABORT("The cutoff values for the multi-grids are not ordered from large to small")
     983        19234 :          ELSE IF (qs_control%e_cutoff(igrid_level) == qs_control%e_cutoff(igrid_level - 1)) THEN
     984            0 :             CPABORT("The same cutoff value was specified for two multi-grids")
     985              :          END IF
     986              :       END DO
     987         9032 :       CALL timestop(handle)
     988        27096 :    END SUBROUTINE read_mgrid_section
     989              : 
     990              : ! **************************************************************************************************
     991              : !> \brief ...
     992              : !> \param qs_control ...
     993              : !> \param qs_section ...
     994              : !> \param cell optional cell used to transform Cartesian input vectors
     995              : ! **************************************************************************************************
     996       144512 :    SUBROUTINE read_qs_section(qs_control, qs_section, cell)
     997              : 
     998              :       TYPE(qs_control_type), INTENT(INOUT)               :: qs_control
     999              :       TYPE(section_vals_type), POINTER                   :: qs_section
    1000              :       TYPE(cell_type), OPTIONAL, POINTER                 :: cell
    1001              : 
    1002              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'read_qs_section'
    1003              : 
    1004              :       CHARACTER(LEN=2)                                   :: element_symbol
    1005              :       CHARACTER(LEN=default_string_length)               :: cval
    1006              :       CHARACTER(LEN=default_string_length), &
    1007         9032 :          DIMENSION(:), POINTER                           :: clist
    1008              :       INTEGER                                            :: handle, itmp, j, jj, k, n_rep, n_var, &
    1009              :                                                             ngauss, ngp, nrep, znum
    1010         9032 :       INTEGER, DIMENSION(:), POINTER                     :: tmplist
    1011              :       LOGICAL :: dftb_scc_mixer_explicit, dftb_tblite_mixer_explicit, explicit, &
    1012              :          tblite_reference_cli, tblite_reference_cli_section, tblite_section_active, was_present, &
    1013              :          xtb_scc_mixer_explicit, xtb_tblite_mixer_explicit
    1014              :       REAL(dp)                                           :: tmp, tmpsqrt, value
    1015         9032 :       REAL(dp), POINTER                                  :: scal(:)
    1016              :       TYPE(section_vals_type), POINTER :: cdft_control_section, ddapc_restraint_section, &
    1017              :          dftb_parameter, dftb_section, dftb_tblite_mixer, eeq_section, genpot_section, &
    1018              :          lri_optbas_section, mull_section, nonbonded_section, s2_restraint_section, se_section, &
    1019              :          xtb_parameter, xtb_section, xtb_tblite, xtb_tblite_mixer, xtb_tblite_ref_cli
    1020              : 
    1021         9032 :       CALL timeset(routineN, handle)
    1022              : 
    1023         9032 :       was_present = .FALSE.
    1024         9032 :       NULLIFY (mull_section, ddapc_restraint_section, s2_restraint_section, &
    1025         9032 :                se_section, dftb_section, xtb_section, dftb_parameter, xtb_parameter, lri_optbas_section, &
    1026         9032 :                cdft_control_section, genpot_section, eeq_section, dftb_tblite_mixer, &
    1027         9032 :                xtb_tblite_mixer, xtb_tblite_ref_cli)
    1028              : 
    1029         9032 :       mull_section => section_vals_get_subs_vals(qs_section, "MULLIKEN_RESTRAINT")
    1030         9032 :       ddapc_restraint_section => section_vals_get_subs_vals(qs_section, "DDAPC_RESTRAINT")
    1031         9032 :       s2_restraint_section => section_vals_get_subs_vals(qs_section, "S2_RESTRAINT")
    1032         9032 :       se_section => section_vals_get_subs_vals(qs_section, "SE")
    1033         9032 :       dftb_section => section_vals_get_subs_vals(qs_section, "DFTB")
    1034         9032 :       xtb_section => section_vals_get_subs_vals(qs_section, "xTB")
    1035         9032 :       dftb_parameter => section_vals_get_subs_vals(dftb_section, "PARAMETER")
    1036         9032 :       dftb_tblite_mixer => section_vals_get_subs_vals(dftb_section, "TBLITE_MIXER")
    1037         9032 :       xtb_parameter => section_vals_get_subs_vals(xtb_section, "PARAMETER")
    1038         9032 :       eeq_section => section_vals_get_subs_vals(xtb_section, "EEQ")
    1039         9032 :       lri_optbas_section => section_vals_get_subs_vals(qs_section, "OPTIMIZE_LRI_BASIS")
    1040         9032 :       cdft_control_section => section_vals_get_subs_vals(qs_section, "CDFT")
    1041         9032 :       nonbonded_section => section_vals_get_subs_vals(xtb_section, "NONBONDED")
    1042         9032 :       genpot_section => section_vals_get_subs_vals(nonbonded_section, "GENPOT")
    1043         9032 :       xtb_tblite_mixer => section_vals_get_subs_vals(xtb_section, "TBLITE_MIXER")
    1044         9032 :       xtb_tblite => section_vals_get_subs_vals(xtb_section, "TBLITE")
    1045         9032 :       xtb_tblite_ref_cli => section_vals_get_subs_vals(xtb_tblite, "REFERENCE_CLI")
    1046              : 
    1047              :       ! Setup all defaults values and overwrite input parameters
    1048              :       ! EPS_DEFAULT should set the target accuracy in the total energy (~per electron) or a closely related value
    1049         9032 :       CALL section_vals_val_get(qs_section, "EPS_DEFAULT", r_val=value)
    1050         9032 :       tmpsqrt = SQRT(value) ! a trick to work around a NAG 5.1 optimizer bug
    1051              : 
    1052              :       ! random choice ?
    1053         9032 :       qs_control%eps_core_charge = value/100.0_dp
    1054              :       ! correct if all Gaussians would have the same radius (overlap will be smaller than eps_pgf_orb**2).
    1055              :       ! Can be significantly in error if not... requires fully new screening/pairlist procedures
    1056         9032 :       qs_control%eps_pgf_orb = tmpsqrt
    1057         9032 :       qs_control%eps_kg_orb = qs_control%eps_pgf_orb
    1058              :       ! consistent since also a kind of overlap
    1059         9032 :       qs_control%eps_ppnl = qs_control%eps_pgf_orb/100.0_dp
    1060              :       ! accuracy is basically set by the overlap, this sets an empirical shift
    1061         9032 :       qs_control%eps_ppl = 1.0E-2_dp
    1062              :       !
    1063         9032 :       qs_control%gapw_control%eps_cpc = value
    1064              :       ! expexted error in the density
    1065         9032 :       qs_control%eps_rho_gspace = value
    1066         9032 :       qs_control%eps_rho_rspace = value
    1067              :       ! error in the gradient, can be the sqrt of the error in the energy, ignored if map_consistent
    1068         9032 :       qs_control%eps_gvg_rspace = tmpsqrt
    1069              :       !
    1070         9032 :       CALL section_vals_val_get(qs_section, "EPS_CORE_CHARGE", n_rep_val=n_rep)
    1071         9032 :       IF (n_rep /= 0) THEN
    1072            0 :          CALL section_vals_val_get(qs_section, "EPS_CORE_CHARGE", r_val=qs_control%eps_core_charge)
    1073              :       END IF
    1074         9032 :       CALL section_vals_val_get(qs_section, "EPS_GVG_RSPACE", n_rep_val=n_rep)
    1075         9032 :       IF (n_rep /= 0) THEN
    1076          164 :          CALL section_vals_val_get(qs_section, "EPS_GVG_RSPACE", r_val=qs_control%eps_gvg_rspace)
    1077              :       END IF
    1078         9032 :       CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", n_rep_val=n_rep)
    1079         9032 :       IF (n_rep /= 0) THEN
    1080          658 :          CALL section_vals_val_get(qs_section, "EPS_PGF_ORB", r_val=qs_control%eps_pgf_orb)
    1081              :       END IF
    1082         9032 :       CALL section_vals_val_get(qs_section, "EPS_KG_ORB", n_rep_val=n_rep)
    1083         9032 :       IF (n_rep /= 0) THEN
    1084           66 :          CALL section_vals_val_get(qs_section, "EPS_KG_ORB", r_val=tmp)
    1085           66 :          qs_control%eps_kg_orb = SQRT(tmp)
    1086              :       END IF
    1087         9032 :       CALL section_vals_val_get(qs_section, "EPS_PPL", n_rep_val=n_rep)
    1088         9032 :       IF (n_rep /= 0) THEN
    1089         9032 :          CALL section_vals_val_get(qs_section, "EPS_PPL", r_val=qs_control%eps_ppl)
    1090              :       END IF
    1091         9032 :       CALL section_vals_val_get(qs_section, "EPS_PPNL", n_rep_val=n_rep)
    1092         9032 :       IF (n_rep /= 0) THEN
    1093            0 :          CALL section_vals_val_get(qs_section, "EPS_PPNL", r_val=qs_control%eps_ppnl)
    1094              :       END IF
    1095         9032 :       CALL section_vals_val_get(qs_section, "EPS_RHO", n_rep_val=n_rep)
    1096         9032 :       IF (n_rep /= 0) THEN
    1097           30 :          CALL section_vals_val_get(qs_section, "EPS_RHO", r_val=qs_control%eps_rho_gspace)
    1098           30 :          qs_control%eps_rho_rspace = qs_control%eps_rho_gspace
    1099              :       END IF
    1100         9032 :       CALL section_vals_val_get(qs_section, "EPS_RHO_RSPACE", n_rep_val=n_rep)
    1101         9032 :       IF (n_rep /= 0) THEN
    1102            2 :          CALL section_vals_val_get(qs_section, "EPS_RHO_RSPACE", r_val=qs_control%eps_rho_rspace)
    1103              :       END IF
    1104         9032 :       CALL section_vals_val_get(qs_section, "EPS_RHO_GSPACE", n_rep_val=n_rep)
    1105         9032 :       IF (n_rep /= 0) THEN
    1106            2 :          CALL section_vals_val_get(qs_section, "EPS_RHO_GSPACE", r_val=qs_control%eps_rho_gspace)
    1107              :       END IF
    1108         9032 :       CALL section_vals_val_get(qs_section, "EPS_FILTER_MATRIX", n_rep_val=n_rep)
    1109         9032 :       IF (n_rep /= 0) THEN
    1110         9032 :          CALL section_vals_val_get(qs_section, "EPS_FILTER_MATRIX", r_val=qs_control%eps_filter_matrix)
    1111              :       END IF
    1112         9032 :       CALL section_vals_val_get(qs_section, "EPS_CPC", n_rep_val=n_rep)
    1113         9032 :       IF (n_rep /= 0) THEN
    1114            0 :          CALL section_vals_val_get(qs_section, "EPS_CPC", r_val=qs_control%gapw_control%eps_cpc)
    1115              :       END IF
    1116              : 
    1117         9032 :       CALL section_vals_val_get(qs_section, "EPSFIT", r_val=qs_control%gapw_control%eps_fit)
    1118         9032 :       CALL section_vals_val_get(qs_section, "EPSISO", r_val=qs_control%gapw_control%eps_iso)
    1119         9032 :       CALL section_vals_val_get(qs_section, "EPSSVD", r_val=qs_control%gapw_control%eps_svd)
    1120         9032 :       CALL section_vals_val_get(qs_section, "EPSRHO0", r_val=qs_control%gapw_control%eps_Vrho0)
    1121         9032 :       CALL section_vals_val_get(qs_section, "ALPHA0_HARD", r_val=qs_control%gapw_control%alpha0_hard)
    1122         9032 :       qs_control%gapw_control%alpha0_hard_from_input = .FALSE.
    1123         9032 :       IF (qs_control%gapw_control%alpha0_hard /= 0.0_dp) qs_control%gapw_control%alpha0_hard_from_input = .TRUE.
    1124         9032 :       CALL section_vals_val_get(qs_section, "FORCE_PAW", l_val=qs_control%gapw_control%force_paw)
    1125         9032 :       CALL section_vals_val_get(qs_section, "MAX_RAD_LOCAL", r_val=qs_control%gapw_control%max_rad_local)
    1126              : 
    1127         9032 :       CALL section_vals_val_get(qs_section, "MIN_PAIR_LIST_RADIUS", r_val=qs_control%pairlist_radius)
    1128              : 
    1129         9032 :       CALL section_vals_val_get(qs_section, "LS_SCF", l_val=qs_control%do_ls_scf)
    1130         9032 :       CALL section_vals_val_get(qs_section, "ALMO_SCF", l_val=qs_control%do_almo_scf)
    1131         9032 :       CALL section_vals_val_get(qs_section, "KG_METHOD", l_val=qs_control%do_kg)
    1132              : 
    1133              :       ! Logicals
    1134         9032 :       CALL section_vals_val_get(qs_section, "REF_EMBED_SUBSYS", l_val=qs_control%ref_embed_subsys)
    1135         9032 :       CALL section_vals_val_get(qs_section, "CLUSTER_EMBED_SUBSYS", l_val=qs_control%cluster_embed_subsys)
    1136         9032 :       CALL section_vals_val_get(qs_section, "HIGH_LEVEL_EMBED_SUBSYS", l_val=qs_control%high_level_embed_subsys)
    1137         9032 :       CALL section_vals_val_get(qs_section, "DFET_EMBEDDED", l_val=qs_control%dfet_embedded)
    1138         9032 :       CALL section_vals_val_get(qs_section, "DMFET_EMBEDDED", l_val=qs_control%dmfet_embedded)
    1139              : 
    1140              :       ! Integers gapw
    1141         9032 :       CALL section_vals_val_get(qs_section, "LMAXN1", i_val=qs_control%gapw_control%lmax_sphere)
    1142         9032 :       CALL section_vals_val_get(qs_section, "LMAXN0", i_val=qs_control%gapw_control%lmax_rho0)
    1143         9032 :       CALL section_vals_val_get(qs_section, "LADDN0", i_val=qs_control%gapw_control%ladd_rho0)
    1144         9032 :       CALL section_vals_val_get(qs_section, "QUADRATURE", i_val=qs_control%gapw_control%quadrature)
    1145              :       ! GAPW 1c basis
    1146         9032 :       CALL section_vals_val_get(qs_section, "GAPW_1C_BASIS", i_val=qs_control%gapw_control%basis_1c)
    1147         9032 :       IF (qs_control%gapw_control%basis_1c /= gapw_1c_orb) THEN
    1148          106 :          qs_control%gapw_control%eps_svd = MAX(qs_control%gapw_control%eps_svd, 1.E-12_dp)
    1149              :       END IF
    1150              :       ! GAPW accurate integration
    1151         9032 :       CALL section_vals_val_get(qs_section, "GAPW_ACCURATE_XCINT", l_val=qs_control%gapw_control%accurate_xcint)
    1152         9032 :       CALL section_vals_val_get(qs_section, "ALPHA_WEIGHTS", r_val=qs_control%gapw_control%aweights)
    1153              : 
    1154              :       ! Integers grids
    1155         9032 :       CALL section_vals_val_get(qs_section, "PW_GRID", i_val=itmp)
    1156            0 :       SELECT CASE (itmp)
    1157              :       CASE (do_pwgrid_spherical)
    1158            0 :          qs_control%pw_grid_opt%spherical = .TRUE.
    1159            0 :          qs_control%pw_grid_opt%fullspace = .FALSE.
    1160              :       CASE (do_pwgrid_ns_fullspace)
    1161         9032 :          qs_control%pw_grid_opt%spherical = .FALSE.
    1162         9032 :          qs_control%pw_grid_opt%fullspace = .TRUE.
    1163              :       CASE (do_pwgrid_ns_halfspace)
    1164            0 :          qs_control%pw_grid_opt%spherical = .FALSE.
    1165         9032 :          qs_control%pw_grid_opt%fullspace = .FALSE.
    1166              :       END SELECT
    1167              : 
    1168              :       !   Method for PPL calculation
    1169         9032 :       CALL section_vals_val_get(qs_section, "CORE_PPL", i_val=itmp)
    1170         9032 :       qs_control%do_ppl_method = itmp
    1171              : 
    1172         9032 :       CALL section_vals_val_get(qs_section, "PW_GRID_LAYOUT", i_vals=tmplist)
    1173        27096 :       qs_control%pw_grid_opt%distribution_layout = tmplist
    1174         9032 :       CALL section_vals_val_get(qs_section, "PW_GRID_BLOCKED", i_val=qs_control%pw_grid_opt%blocked)
    1175              : 
    1176              :       !Integers extrapolation
    1177         9032 :       CALL section_vals_val_get(qs_section, "EXTRAPOLATION", i_val=qs_control%wf_interpolation_method_nr)
    1178         9032 :       CALL section_vals_val_get(qs_section, "EXTRAPOLATION_ORDER", i_val=qs_control%wf_extrapolation_order)
    1179              : 
    1180              :       !Method
    1181         9032 :       CALL section_vals_val_get(qs_section, "METHOD", i_val=qs_control%method_id)
    1182         9032 :       qs_control%gapw = .FALSE.
    1183         9032 :       qs_control%gapw_xc = .FALSE.
    1184         9032 :       qs_control%gpw = .FALSE.
    1185         9032 :       qs_control%pao = .FALSE.
    1186         9032 :       qs_control%dftb = .FALSE.
    1187         9032 :       qs_control%xtb = .FALSE.
    1188         9032 :       qs_control%semi_empirical = .FALSE.
    1189         9032 :       qs_control%ofgpw = .FALSE.
    1190         9032 :       qs_control%lrigpw = .FALSE.
    1191         9032 :       qs_control%rigpw = .FALSE.
    1192        10306 :       SELECT CASE (qs_control%method_id)
    1193              :       CASE (do_method_gapw)
    1194         1274 :          CALL cite_reference(Lippert1999)
    1195         1274 :          CALL cite_reference(Krack2000)
    1196         1274 :          qs_control%gapw = .TRUE.
    1197              :       CASE (do_method_gapw_xc)
    1198          184 :          qs_control%gapw_xc = .TRUE.
    1199              :       CASE (do_method_gpw)
    1200         5004 :          CALL cite_reference(Lippert1997)
    1201         5004 :          CALL cite_reference(VandeVondele2005a)
    1202         5004 :          qs_control%gpw = .TRUE.
    1203              :       CASE (do_method_ofgpw)
    1204            0 :          qs_control%ofgpw = .TRUE.
    1205              :       CASE (do_method_lrigpw)
    1206           42 :          qs_control%lrigpw = .TRUE.
    1207              :       CASE (do_method_rigpw)
    1208            2 :          qs_control%rigpw = .TRUE.
    1209              :       CASE (do_method_dftb)
    1210          298 :          qs_control%dftb = .TRUE.
    1211          298 :          CALL cite_reference(Porezag1995)
    1212          298 :          CALL cite_reference(Seifert1996)
    1213              :       CASE (do_method_xtb)
    1214         1228 :          qs_control%xtb = .TRUE.
    1215         1228 :          CALL cite_reference(Grimme2017)
    1216         1228 :          CALL cite_reference(Pracht2019)
    1217              :       CASE (do_method_mndo)
    1218           52 :          CALL cite_reference(Dewar1977)
    1219           52 :          qs_control%semi_empirical = .TRUE.
    1220              :       CASE (do_method_am1)
    1221          112 :          CALL cite_reference(Dewar1985)
    1222          112 :          qs_control%semi_empirical = .TRUE.
    1223              :       CASE (do_method_pm3)
    1224           48 :          CALL cite_reference(Stewart1989)
    1225           48 :          qs_control%semi_empirical = .TRUE.
    1226              :       CASE (do_method_pnnl)
    1227           14 :          CALL cite_reference(Schenter2008)
    1228           14 :          qs_control%semi_empirical = .TRUE.
    1229              :       CASE (do_method_pm6)
    1230          754 :          CALL cite_reference(Stewart2007)
    1231          754 :          qs_control%semi_empirical = .TRUE.
    1232              :       CASE (do_method_pm6fm)
    1233            0 :          CALL cite_reference(VanVoorhis2015)
    1234            0 :          qs_control%semi_empirical = .TRUE.
    1235              :       CASE (do_method_pdg)
    1236            2 :          CALL cite_reference(Repasky2002)
    1237            2 :          qs_control%semi_empirical = .TRUE.
    1238              :       CASE (do_method_rm1)
    1239            2 :          CALL cite_reference(Rocha2006)
    1240            2 :          qs_control%semi_empirical = .TRUE.
    1241              :       CASE (do_method_mndod)
    1242           16 :          CALL cite_reference(Dewar1977)
    1243           16 :          CALL cite_reference(Thiel1992)
    1244         9048 :          qs_control%semi_empirical = .TRUE.
    1245              :       END SELECT
    1246              : 
    1247         9032 :       CALL section_vals_get(mull_section, explicit=qs_control%mulliken_restraint)
    1248              : 
    1249         9032 :       IF (qs_control%mulliken_restraint) THEN
    1250            2 :          CALL section_vals_val_get(mull_section, "STRENGTH", r_val=qs_control%mulliken_restraint_control%strength)
    1251            2 :          CALL section_vals_val_get(mull_section, "TARGET", r_val=qs_control%mulliken_restraint_control%target)
    1252            2 :          CALL section_vals_val_get(mull_section, "ATOMS", n_rep_val=n_rep)
    1253            2 :          jj = 0
    1254            4 :          DO k = 1, n_rep
    1255            2 :             CALL section_vals_val_get(mull_section, "ATOMS", i_rep_val=k, i_vals=tmplist)
    1256            4 :             jj = jj + SIZE(tmplist)
    1257              :          END DO
    1258            2 :          qs_control%mulliken_restraint_control%natoms = jj
    1259            2 :          IF (qs_control%mulliken_restraint_control%natoms < 1) THEN
    1260            0 :             CPABORT("Need at least 1 atom to use mulliken constraints")
    1261              :          END IF
    1262            6 :          ALLOCATE (qs_control%mulliken_restraint_control%atoms(qs_control%mulliken_restraint_control%natoms))
    1263            2 :          jj = 0
    1264            6 :          DO k = 1, n_rep
    1265            2 :             CALL section_vals_val_get(mull_section, "ATOMS", i_rep_val=k, i_vals=tmplist)
    1266            6 :             DO j = 1, SIZE(tmplist)
    1267            2 :                jj = jj + 1
    1268            4 :                qs_control%mulliken_restraint_control%atoms(jj) = tmplist(j)
    1269              :             END DO
    1270              :          END DO
    1271              :       END IF
    1272         9032 :       CALL section_vals_get(ddapc_restraint_section, n_repetition=nrep, explicit=qs_control%ddapc_restraint)
    1273         9032 :       IF (qs_control%ddapc_restraint) THEN
    1274           60 :          ALLOCATE (qs_control%ddapc_restraint_control(nrep))
    1275           14 :          CALL read_ddapc_section(qs_control, qs_section=qs_section)
    1276           14 :          qs_control%ddapc_restraint_is_spin = .FALSE.
    1277           14 :          qs_control%ddapc_explicit_potential = .FALSE.
    1278              :       END IF
    1279              : 
    1280         9032 :       CALL section_vals_get(s2_restraint_section, explicit=qs_control%s2_restraint)
    1281         9032 :       IF (qs_control%s2_restraint) THEN
    1282              :          CALL section_vals_val_get(s2_restraint_section, "STRENGTH", &
    1283            0 :                                    r_val=qs_control%s2_restraint_control%strength)
    1284              :          CALL section_vals_val_get(s2_restraint_section, "TARGET", &
    1285            0 :                                    r_val=qs_control%s2_restraint_control%target)
    1286              :          CALL section_vals_val_get(s2_restraint_section, "FUNCTIONAL_FORM", &
    1287            0 :                                    i_val=qs_control%s2_restraint_control%functional_form)
    1288              :       END IF
    1289              : 
    1290         9032 :       CALL section_vals_get(cdft_control_section, explicit=qs_control%cdft)
    1291         9032 :       IF (qs_control%cdft) THEN
    1292          296 :          CALL read_cdft_control_section(qs_control, cdft_control_section)
    1293              :       END IF
    1294              : 
    1295              :       ! Semi-empirical code
    1296         9032 :       IF (qs_control%semi_empirical) THEN
    1297              :          CALL section_vals_val_get(se_section, "ORTHOGONAL_BASIS", &
    1298         1000 :                                    l_val=qs_control%se_control%orthogonal_basis)
    1299              :          CALL section_vals_val_get(se_section, "DELTA", &
    1300         1000 :                                    r_val=qs_control%se_control%delta)
    1301              :          CALL section_vals_val_get(se_section, "ANALYTICAL_GRADIENTS", &
    1302         1000 :                                    l_val=qs_control%se_control%analytical_gradients)
    1303              :          CALL section_vals_val_get(se_section, "FORCE_KDSO-D_EXCHANGE", &
    1304         1000 :                                    l_val=qs_control%se_control%force_kdsod_EX)
    1305              :          ! Integral Screening
    1306              :          CALL section_vals_val_get(se_section, "INTEGRAL_SCREENING", &
    1307         1000 :                                    i_val=qs_control%se_control%integral_screening)
    1308         1000 :          IF (qs_control%method_id == do_method_pnnl) THEN
    1309           14 :             IF (qs_control%se_control%integral_screening /= do_se_IS_slater) THEN
    1310              :                CALL cp_warn(__LOCATION__, &
    1311              :                             "PNNL semi-empirical parameterization supports only the Slater type "// &
    1312            0 :                             "integral scheme. Revert to Slater and continue the calculation.")
    1313              :             END IF
    1314           14 :             qs_control%se_control%integral_screening = do_se_IS_slater
    1315              :          END IF
    1316              :          ! Global Arrays variable
    1317              :          CALL section_vals_val_get(se_section, "GA%NCELLS", &
    1318         1000 :                                    i_val=qs_control%se_control%ga_ncells)
    1319              :          ! Long-Range correction
    1320              :          CALL section_vals_val_get(se_section, "LR_CORRECTION%CUTOFF", &
    1321         1000 :                                    r_val=qs_control%se_control%cutoff_lrc)
    1322         1000 :          qs_control%se_control%taper_lrc = qs_control%se_control%cutoff_lrc
    1323              :          CALL section_vals_val_get(se_section, "LR_CORRECTION%RC_TAPER", &
    1324         1000 :                                    explicit=explicit)
    1325         1000 :          IF (explicit) THEN
    1326              :             CALL section_vals_val_get(se_section, "LR_CORRECTION%RC_TAPER", &
    1327            0 :                                       r_val=qs_control%se_control%taper_lrc)
    1328              :          END IF
    1329              :          CALL section_vals_val_get(se_section, "LR_CORRECTION%RC_RANGE", &
    1330         1000 :                                    r_val=qs_control%se_control%range_lrc)
    1331              :          ! Coulomb
    1332              :          CALL section_vals_val_get(se_section, "COULOMB%CUTOFF", &
    1333         1000 :                                    r_val=qs_control%se_control%cutoff_cou)
    1334         1000 :          qs_control%se_control%taper_cou = qs_control%se_control%cutoff_cou
    1335              :          CALL section_vals_val_get(se_section, "COULOMB%RC_TAPER", &
    1336         1000 :                                    explicit=explicit)
    1337         1000 :          IF (explicit) THEN
    1338              :             CALL section_vals_val_get(se_section, "COULOMB%RC_TAPER", &
    1339            0 :                                       r_val=qs_control%se_control%taper_cou)
    1340              :          END IF
    1341              :          CALL section_vals_val_get(se_section, "COULOMB%RC_RANGE", &
    1342         1000 :                                    r_val=qs_control%se_control%range_cou)
    1343              :          ! Exchange
    1344              :          CALL section_vals_val_get(se_section, "EXCHANGE%CUTOFF", &
    1345         1000 :                                    r_val=qs_control%se_control%cutoff_exc)
    1346         1000 :          qs_control%se_control%taper_exc = qs_control%se_control%cutoff_exc
    1347              :          CALL section_vals_val_get(se_section, "EXCHANGE%RC_TAPER", &
    1348         1000 :                                    explicit=explicit)
    1349         1000 :          IF (explicit) THEN
    1350              :             CALL section_vals_val_get(se_section, "EXCHANGE%RC_TAPER", &
    1351           38 :                                       r_val=qs_control%se_control%taper_exc)
    1352              :          END IF
    1353              :          CALL section_vals_val_get(se_section, "EXCHANGE%RC_RANGE", &
    1354         1000 :                                    r_val=qs_control%se_control%range_exc)
    1355              :          ! Screening (only if the integral scheme is of dumped type)
    1356         1000 :          IF (qs_control%se_control%integral_screening == do_se_IS_kdso_d) THEN
    1357              :             CALL section_vals_val_get(se_section, "SCREENING%RC_TAPER", &
    1358           14 :                                       r_val=qs_control%se_control%taper_scr)
    1359              :             CALL section_vals_val_get(se_section, "SCREENING%RC_RANGE", &
    1360           14 :                                       r_val=qs_control%se_control%range_scr)
    1361              :          END IF
    1362              :          ! Periodic Type Calculation
    1363              :          CALL section_vals_val_get(se_section, "PERIODIC", &
    1364         1000 :                                    i_val=qs_control%se_control%periodic_type)
    1365         1968 :          SELECT CASE (qs_control%se_control%periodic_type)
    1366              :          CASE (do_se_lr_none)
    1367          968 :             qs_control%se_control%do_ewald = .FALSE.
    1368          968 :             qs_control%se_control%do_ewald_r3 = .FALSE.
    1369          968 :             qs_control%se_control%do_ewald_gks = .FALSE.
    1370              :          CASE (do_se_lr_ewald)
    1371           30 :             qs_control%se_control%do_ewald = .TRUE.
    1372           30 :             qs_control%se_control%do_ewald_r3 = .FALSE.
    1373           30 :             qs_control%se_control%do_ewald_gks = .FALSE.
    1374              :          CASE (do_se_lr_ewald_gks)
    1375            2 :             qs_control%se_control%do_ewald = .FALSE.
    1376            2 :             qs_control%se_control%do_ewald_r3 = .FALSE.
    1377            2 :             qs_control%se_control%do_ewald_gks = .TRUE.
    1378            2 :             IF (qs_control%method_id /= do_method_pnnl) THEN
    1379              :                CALL cp_abort(__LOCATION__, &
    1380              :                              "A periodic semi-empirical calculation was requested with a long-range  "// &
    1381              :                              "summation on the single integral evaluation. This scheme is supported  "// &
    1382            0 :                              "only by the PNNL parameterization.")
    1383              :             END IF
    1384              :          CASE (do_se_lr_ewald_r3)
    1385            0 :             qs_control%se_control%do_ewald = .TRUE.
    1386            0 :             qs_control%se_control%do_ewald_r3 = .TRUE.
    1387            0 :             qs_control%se_control%do_ewald_gks = .FALSE.
    1388         1000 :             IF (qs_control%se_control%integral_screening /= do_se_IS_kdso) THEN
    1389              :                CALL cp_abort(__LOCATION__, &
    1390              :                              "A periodic semi-empirical calculation was requested with a long-range  "// &
    1391              :                              "summation for the slowly convergent part 1/R^3, which is not congruent "// &
    1392              :                              "with the integral screening chosen. The only integral screening supported "// &
    1393            0 :                              "by this periodic type calculation is the standard Klopman-Dewar-Sabelli-Ohno.")
    1394              :             END IF
    1395              :          END SELECT
    1396              : 
    1397              :          ! dispersion pair potentials
    1398              :          CALL section_vals_val_get(se_section, "DISPERSION", &
    1399         1000 :                                    l_val=qs_control%se_control%dispersion)
    1400              :          CALL section_vals_val_get(se_section, "DISPERSION_RADIUS", &
    1401         1000 :                                    r_val=qs_control%se_control%rcdisp)
    1402              :          CALL section_vals_val_get(se_section, "COORDINATION_CUTOFF", &
    1403         1000 :                                    r_val=qs_control%se_control%epscn)
    1404         1000 :          CALL section_vals_val_get(se_section, "D3_SCALING", r_vals=scal)
    1405         1000 :          qs_control%se_control%sd3(1) = scal(1)
    1406         1000 :          qs_control%se_control%sd3(2) = scal(2)
    1407         1000 :          qs_control%se_control%sd3(3) = scal(3)
    1408              :          CALL section_vals_val_get(se_section, "DISPERSION_PARAMETER_FILE", &
    1409         1000 :                                    c_val=qs_control%se_control%dispersion_parameter_file)
    1410              : 
    1411              :          ! Stop the execution for non-implemented features
    1412         1000 :          IF (qs_control%se_control%periodic_type == do_se_lr_ewald_r3) THEN
    1413            0 :             CPABORT("EWALD_R3 not implemented yet!")
    1414              :          END IF
    1415              : 
    1416              :          IF (qs_control%method_id == do_method_mndo .OR. &
    1417              :              qs_control%method_id == do_method_am1 .OR. &
    1418              :              qs_control%method_id == do_method_mndod .OR. &
    1419              :              qs_control%method_id == do_method_pdg .OR. &
    1420              :              qs_control%method_id == do_method_pm3 .OR. &
    1421              :              qs_control%method_id == do_method_pm6 .OR. &
    1422              :              qs_control%method_id == do_method_pm6fm .OR. &
    1423         1000 :              qs_control%method_id == do_method_pnnl .OR. &
    1424              :              qs_control%method_id == do_method_rm1) THEN
    1425         1000 :             qs_control%se_control%orthogonal_basis = .TRUE.
    1426              :          END IF
    1427              :       END IF
    1428              : 
    1429              :       ! DFTB code
    1430         9032 :       IF (qs_control%dftb) THEN
    1431              :          CALL section_vals_val_get(dftb_section, "ORTHOGONAL_BASIS", &
    1432          298 :                                    l_val=qs_control%dftb_control%orthogonal_basis)
    1433              :          CALL section_vals_val_get(dftb_section, "SELF_CONSISTENT", &
    1434          298 :                                    l_val=qs_control%dftb_control%self_consistent)
    1435              :          CALL section_vals_val_get(dftb_section, "DISPERSION", &
    1436          298 :                                    l_val=qs_control%dftb_control%dispersion)
    1437              :          CALL section_vals_val_get(dftb_section, "DIAGONAL_DFTB3", &
    1438          298 :                                    l_val=qs_control%dftb_control%dftb3_diagonal)
    1439              :          CALL section_vals_val_get(dftb_section, "HB_SR_GAMMA", &
    1440          298 :                                    l_val=qs_control%dftb_control%hb_sr_damp)
    1441              :          CALL section_vals_val_get(dftb_section, "SCC_MIXER", &
    1442          298 :                                    explicit=dftb_scc_mixer_explicit)
    1443              :          CALL section_vals_val_get(dftb_section, "SCC_MIXER", &
    1444          298 :                                    i_val=qs_control%dftb_control%tblite_scc_mixer)
    1445          298 :          CALL section_vals_get(dftb_tblite_mixer, explicit=dftb_tblite_mixer_explicit)
    1446              :          CALL read_tblite_mixer_section(dftb_tblite_mixer, &
    1447              :                                         qs_control%dftb_control%tblite_mixer_iterations, &
    1448              :                                         qs_control%dftb_control%tblite_mixer_memory, &
    1449              :                                         qs_control%dftb_control%tblite_mixer_solver, &
    1450              :                                         qs_control%dftb_control%tblite_mixer_damping, &
    1451              :                                         qs_control%dftb_control%tblite_mixer_omega0, &
    1452              :                                         qs_control%dftb_control%tblite_mixer_min_weight, &
    1453              :                                         qs_control%dftb_control%tblite_mixer_max_weight, &
    1454              :                                         qs_control%dftb_control%tblite_mixer_weight_factor, &
    1455          298 :                                         "DFTB/TBLITE_MIXER")
    1456          298 :          IF (qs_control%do_ls_scf) THEN
    1457           44 :             IF (dftb_scc_mixer_explicit .AND. &
    1458              :                 qs_control%dftb_control%tblite_scc_mixer /= tblite_scc_mixer_none) THEN
    1459              :                CALL cp_warn(__LOCATION__, &
    1460              :                             "DFTB/SCC_MIXER is reset to NONE with QS/LS_SCF; LS_SCF optimizes "// &
    1461            2 :                             "the density matrix directly.")
    1462              :             END IF
    1463           44 :             IF (dftb_tblite_mixer_explicit) THEN
    1464              :                CALL cp_warn(__LOCATION__, &
    1465              :                             "DFTB/TBLITE_MIXER settings are ignored with QS/LS_SCF; LS_SCF controls "// &
    1466            0 :                             "the density-matrix optimization.")
    1467              :             END IF
    1468           44 :             qs_control%dftb_control%tblite_scc_mixer = tblite_scc_mixer_none
    1469              :          END IF
    1470          298 :          IF (qs_control%dftb_control%tblite_mixer_damping <= 0.0_dp) THEN
    1471            0 :             CPABORT("DFTB/TBLITE_MIXER/DAMPING must be positive")
    1472              :          END IF
    1473              :          CALL section_vals_val_get(dftb_section, "EPS_DISP", &
    1474          298 :                                    r_val=qs_control%dftb_control%eps_disp)
    1475          298 :          CALL section_vals_val_get(dftb_section, "DO_EWALD", explicit=explicit)
    1476          298 :          IF (explicit) THEN
    1477              :             CALL section_vals_val_get(dftb_section, "DO_EWALD", &
    1478          206 :                                       l_val=qs_control%dftb_control%do_ewald)
    1479              :          ELSE
    1480           92 :             qs_control%dftb_control%do_ewald = (qs_control%periodicity /= 0)
    1481              :          END IF
    1482              :          CALL section_vals_val_get(dftb_parameter, "PARAM_FILE_PATH", &
    1483          298 :                                    c_val=qs_control%dftb_control%sk_file_path)
    1484              :          CALL section_vals_val_get(dftb_parameter, "PARAM_FILE_NAME", &
    1485          298 :                                    c_val=qs_control%dftb_control%sk_file_list)
    1486              :          CALL section_vals_val_get(dftb_parameter, "HB_SR_PARAM", &
    1487          298 :                                    r_val=qs_control%dftb_control%hb_sr_para)
    1488          298 :          CALL section_vals_val_get(dftb_parameter, "SK_FILE", n_rep_val=n_var)
    1489          644 :          ALLOCATE (qs_control%dftb_control%sk_pair_list(3, n_var))
    1490          394 :          DO k = 1, n_var
    1491              :             CALL section_vals_val_get(dftb_parameter, "SK_FILE", i_rep_val=k, &
    1492           96 :                                       c_vals=clist)
    1493          682 :             qs_control%dftb_control%sk_pair_list(1:3, k) = clist(1:3)
    1494              :          END DO
    1495              :          ! Dispersion type
    1496              :          CALL section_vals_val_get(dftb_parameter, "DISPERSION_TYPE", &
    1497          298 :                                    i_val=qs_control%dftb_control%dispersion_type)
    1498              :          CALL section_vals_val_get(dftb_parameter, "UFF_FORCE_FIELD", &
    1499          298 :                                    c_val=qs_control%dftb_control%uff_force_field)
    1500              :          ! D3 Dispersion
    1501              :          CALL section_vals_val_get(dftb_parameter, "DISPERSION_RADIUS", &
    1502          298 :                                    r_val=qs_control%dftb_control%rcdisp)
    1503              :          CALL section_vals_val_get(dftb_parameter, "COORDINATION_CUTOFF", &
    1504          298 :                                    r_val=qs_control%dftb_control%epscn)
    1505              :          CALL section_vals_val_get(dftb_parameter, "D2_EXP_PRE", &
    1506          298 :                                    r_val=qs_control%dftb_control%exp_pre)
    1507              :          CALL section_vals_val_get(dftb_parameter, "D2_SCALING", &
    1508          298 :                                    r_val=qs_control%dftb_control%scaling)
    1509          298 :          CALL section_vals_val_get(dftb_parameter, "D3_SCALING", r_vals=scal)
    1510          298 :          qs_control%dftb_control%sd3(1) = scal(1)
    1511          298 :          qs_control%dftb_control%sd3(2) = scal(2)
    1512          298 :          qs_control%dftb_control%sd3(3) = scal(3)
    1513          298 :          CALL section_vals_val_get(dftb_parameter, "D3BJ_SCALING", r_vals=scal)
    1514          298 :          qs_control%dftb_control%sd3bj(1) = scal(1)
    1515          298 :          qs_control%dftb_control%sd3bj(2) = scal(2)
    1516          298 :          qs_control%dftb_control%sd3bj(3) = scal(3)
    1517          298 :          qs_control%dftb_control%sd3bj(4) = scal(4)
    1518              :          CALL section_vals_val_get(dftb_parameter, "DISPERSION_PARAMETER_FILE", &
    1519          298 :                                    c_val=qs_control%dftb_control%dispersion_parameter_file)
    1520              : 
    1521          298 :          IF (qs_control%dftb_control%dispersion) CALL cite_reference(Zhechkov2005)
    1522          298 :          IF (qs_control%dftb_control%self_consistent) CALL cite_reference(Elstner1998)
    1523         1490 :          IF (qs_control%dftb_control%hb_sr_damp) CALL cite_reference(Hu2007)
    1524              :       END IF
    1525              : 
    1526              :       ! xTB code
    1527         9032 :       IF (qs_control%xtb) THEN
    1528         1228 :          CALL section_vals_val_get(xtb_section, "GFN_TYPE", i_val=qs_control%xtb_control%gfn_type)
    1529         1228 :          CALL section_vals_val_get(xtb_tblite, "_SECTION_PARAMETERS_", l_val=tblite_section_active)
    1530         1228 :          qs_control%xtb_control%do_tblite = (qs_control%xtb_control%gfn_type == gfn_tblite)
    1531         1228 :          IF (qs_control%xtb_control%do_tblite) THEN
    1532          188 :             IF (.NOT. tblite_section_active) THEN
    1533            0 :                CPABORT("XTB/GFN_TYPE TBLITE requires an XTB/TBLITE section")
    1534              :             END IF
    1535              :             ! The CP2K-internal GFN1 defaults are still used to initialize shared xTB fields.
    1536          188 :             qs_control%xtb_control%gfn_type = gfn1xtb
    1537         1040 :          ELSE IF (tblite_section_active) THEN
    1538            0 :             CPABORT("The XTB/TBLITE section requires XTB/GFN_TYPE TBLITE")
    1539              :          END IF
    1540              :          CALL section_vals_val_get(xtb_section, "SCC_MIXER", &
    1541         1228 :                                    explicit=xtb_scc_mixer_explicit)
    1542              :          CALL section_vals_val_get(xtb_section, "SCC_MIXER", &
    1543         1228 :                                    i_val=qs_control%xtb_control%tblite_scc_mixer)
    1544         1228 :          CALL section_vals_get(xtb_tblite_mixer, explicit=xtb_tblite_mixer_explicit)
    1545              :          CALL read_tblite_mixer_section(xtb_tblite_mixer, &
    1546              :                                         qs_control%xtb_control%tblite_mixer_iterations, &
    1547              :                                         qs_control%xtb_control%tblite_mixer_memory, &
    1548              :                                         qs_control%xtb_control%tblite_mixer_solver, &
    1549              :                                         qs_control%xtb_control%tblite_mixer_damping, &
    1550              :                                         qs_control%xtb_control%tblite_mixer_omega0, &
    1551              :                                         qs_control%xtb_control%tblite_mixer_min_weight, &
    1552              :                                         qs_control%xtb_control%tblite_mixer_max_weight, &
    1553              :                                         qs_control%xtb_control%tblite_mixer_weight_factor, &
    1554         1228 :                                         "XTB/TBLITE_MIXER")
    1555         1228 :          IF (xtb_tblite_mixer_explicit) THEN
    1556              :             CALL section_vals_val_get(xtb_tblite_mixer, "DAMPING", &
    1557            2 :                                       explicit=qs_control%xtb_control%tblite_mixer_damping_explicit)
    1558              :          END IF
    1559         1228 :          IF ((.NOT. qs_control%xtb_control%do_tblite) .AND. &
    1560              :              qs_control%xtb_control%gfn_type == 0) THEN
    1561              :             IF (xtb_scc_mixer_explicit .AND. &
    1562          694 :                 qs_control%xtb_control%tblite_scc_mixer /= tblite_scc_mixer_auto .AND. &
    1563              :                 qs_control%xtb_control%tblite_scc_mixer /= tblite_scc_mixer_none) THEN
    1564              :                CALL cp_warn(__LOCATION__, &
    1565              :                             "XTB/SCC_MIXER is reset to NONE for CP2K-internal GFN0-xTB; "// &
    1566            0 :                             "GFN0-xTB has no SCC variables to mix.")
    1567              :             END IF
    1568          694 :             IF (xtb_tblite_mixer_explicit) THEN
    1569              :                CALL cp_warn(__LOCATION__, &
    1570              :                             "XTB/TBLITE_MIXER settings are ignored for CP2K-internal GFN0-xTB; "// &
    1571            0 :                             "GFN0-xTB has no SCC variables to mix.")
    1572              :             END IF
    1573          694 :             qs_control%xtb_control%tblite_scc_mixer = tblite_scc_mixer_none
    1574              :          END IF
    1575         1228 :          IF (qs_control%do_ls_scf) THEN
    1576           36 :             IF (xtb_scc_mixer_explicit .AND. &
    1577              :                 qs_control%xtb_control%tblite_scc_mixer /= tblite_scc_mixer_none) THEN
    1578              :                CALL cp_warn(__LOCATION__, &
    1579              :                             "XTB/SCC_MIXER is reset to NONE with QS/LS_SCF; LS_SCF optimizes "// &
    1580            4 :                             "the density matrix directly.")
    1581              :             END IF
    1582           36 :             IF (xtb_tblite_mixer_explicit) THEN
    1583              :                CALL cp_warn(__LOCATION__, &
    1584              :                             "XTB/TBLITE_MIXER settings are ignored with QS/LS_SCF; LS_SCF controls "// &
    1585            0 :                             "the density-matrix optimization.")
    1586              :             END IF
    1587           36 :             qs_control%xtb_control%tblite_scc_mixer = tblite_scc_mixer_none
    1588              :          END IF
    1589         1228 :          IF (qs_control%xtb_control%tblite_mixer_damping <= 0.0_dp) THEN
    1590            0 :             CPABORT("XTB/TBLITE_MIXER/DAMPING must be positive")
    1591              :          END IF
    1592         1228 :          CALL section_vals_val_get(xtb_section, "DO_EWALD", explicit=explicit)
    1593         1228 :          IF (explicit) THEN
    1594              :             CALL section_vals_val_get(xtb_section, "DO_EWALD", &
    1595          776 :                                       l_val=qs_control%xtb_control%do_ewald)
    1596              :          ELSE
    1597          452 :             qs_control%xtb_control%do_ewald = (qs_control%periodicity /= 0)
    1598              :          END IF
    1599              :          ! Spin Polarisation
    1600              :          CALL section_vals_val_get(xtb_section, "SPIN_POLARISATION", &
    1601         1228 :                                    l_val=qs_control%xtb_control%do_spinpol)
    1602              :          ! vdW
    1603         1228 :          CALL section_vals_val_get(xtb_section, "VDW_POTENTIAL", explicit=explicit)
    1604         1228 :          IF (explicit) THEN
    1605          684 :             CALL section_vals_val_get(xtb_section, "VDW_POTENTIAL", c_val=cval)
    1606          684 :             CALL uppercase(cval)
    1607            2 :             SELECT CASE (cval)
    1608              :             CASE ("NONE")
    1609            2 :                qs_control%xtb_control%vdw_type = xtb_vdw_type_none
    1610              :             CASE ("DFTD3")
    1611           56 :                qs_control%xtb_control%vdw_type = xtb_vdw_type_d3
    1612              :             CASE ("DFTD4")
    1613          626 :                qs_control%xtb_control%vdw_type = xtb_vdw_type_d4
    1614              :             CASE DEFAULT
    1615          684 :                CPABORT("vdW type")
    1616              :             END SELECT
    1617              :          ELSE
    1618          560 :             SELECT CASE (qs_control%xtb_control%gfn_type)
    1619              :             CASE (0)
    1620           16 :                qs_control%xtb_control%vdw_type = xtb_vdw_type_d4
    1621              :             CASE (1)
    1622          528 :                qs_control%xtb_control%vdw_type = xtb_vdw_type_d3
    1623              :             CASE (2)
    1624            0 :                qs_control%xtb_control%vdw_type = xtb_vdw_type_d4
    1625            0 :                CPABORT("gfn2-xtb tbd")
    1626              :             CASE DEFAULT
    1627          544 :                CPABORT("GFN type")
    1628              :             END SELECT
    1629              :          END IF
    1630              :          !
    1631         1228 :          CALL section_vals_val_get(xtb_section, "STO_NG", i_val=ngauss)
    1632         1228 :          qs_control%xtb_control%sto_ng = ngauss
    1633         1228 :          CALL section_vals_val_get(xtb_section, "HYDROGEN_STO_NG", i_val=ngauss)
    1634         1228 :          qs_control%xtb_control%h_sto_ng = ngauss
    1635         1228 :          CALL section_vals_val_get(xtb_section, "STO_FLEX", explicit=explicit)
    1636         1228 :          IF (explicit) THEN
    1637              :             CALL section_vals_val_get(xtb_section, "STO_FLEX", &
    1638            6 :                                       l_val=qs_control%xtb_control%sto_flex)
    1639              :          ELSE
    1640         1222 :             qs_control%xtb_control%sto_flex = .FALSE.
    1641              :          END IF
    1642              :          CALL section_vals_val_get(xtb_parameter, "PARAM_FILE_PATH", &
    1643         1228 :                                    c_val=qs_control%xtb_control%parameter_file_path)
    1644         1228 :          CALL section_vals_val_get(xtb_parameter, "PARAM_FILE_NAME", explicit=explicit)
    1645         1228 :          IF (explicit) THEN
    1646              :             CALL section_vals_val_get(xtb_parameter, "PARAM_FILE_NAME", &
    1647            0 :                                       c_val=qs_control%xtb_control%parameter_file_name)
    1648              :          ELSE
    1649         1922 :             SELECT CASE (qs_control%xtb_control%gfn_type)
    1650              :             CASE (0)
    1651          694 :                qs_control%xtb_control%parameter_file_name = "xTB0_parameters"
    1652              :             CASE (1)
    1653          534 :                qs_control%xtb_control%parameter_file_name = "xTB1_parameters"
    1654              :             CASE (2)
    1655            0 :                CPABORT("gfn2-xtb tbd")
    1656              :             CASE DEFAULT
    1657         1228 :                CPABORT("GFN type")
    1658              :             END SELECT
    1659              :          END IF
    1660              :          !
    1661              :          CALL section_vals_val_get(xtb_parameter, "SPINPOL_PARAM_FILE_NAME", &
    1662         1228 :                                    c_val=qs_control%xtb_control%spinpol_param_file_name)
    1663              :          ! D3 Dispersion
    1664              :          CALL section_vals_val_get(xtb_parameter, "DISPERSION_RADIUS", &
    1665         1228 :                                    r_val=qs_control%xtb_control%rcdisp)
    1666              :          CALL section_vals_val_get(xtb_parameter, "COORDINATION_CUTOFF", &
    1667         1228 :                                    r_val=qs_control%xtb_control%epscn)
    1668         1228 :          CALL section_vals_val_get(xtb_parameter, "D3BJ_SCALING", explicit=explicit)
    1669         1228 :          IF (explicit) THEN
    1670            0 :             CALL section_vals_val_get(xtb_parameter, "D3BJ_SCALING", r_vals=scal)
    1671            0 :             qs_control%xtb_control%s6 = scal(1)
    1672            0 :             qs_control%xtb_control%s8 = scal(2)
    1673              :          ELSE
    1674         1922 :             SELECT CASE (qs_control%xtb_control%gfn_type)
    1675              :             CASE (0)
    1676          694 :                qs_control%xtb_control%s6 = 1.00_dp
    1677          694 :                qs_control%xtb_control%s8 = 2.85_dp
    1678              :             CASE (1)
    1679          534 :                qs_control%xtb_control%s6 = 1.00_dp
    1680          534 :                qs_control%xtb_control%s8 = 2.40_dp
    1681              :             CASE (2)
    1682            0 :                CPABORT("gfn2-xtb tbd")
    1683              :             CASE DEFAULT
    1684         1228 :                CPABORT("GFN type")
    1685              :             END SELECT
    1686              :          END IF
    1687         1228 :          CALL section_vals_val_get(xtb_parameter, "D3BJ_PARAM", explicit=explicit)
    1688         1228 :          IF (explicit) THEN
    1689            0 :             CALL section_vals_val_get(xtb_parameter, "D3BJ_PARAM", r_vals=scal)
    1690            0 :             qs_control%xtb_control%a1 = scal(1)
    1691            0 :             qs_control%xtb_control%a2 = scal(2)
    1692              :          ELSE
    1693         1922 :             SELECT CASE (qs_control%xtb_control%gfn_type)
    1694              :             CASE (0)
    1695          694 :                qs_control%xtb_control%a1 = 0.80_dp
    1696          694 :                qs_control%xtb_control%a2 = 4.60_dp
    1697              :             CASE (1)
    1698          534 :                qs_control%xtb_control%a1 = 0.63_dp
    1699          534 :                qs_control%xtb_control%a2 = 5.00_dp
    1700              :             CASE (2)
    1701            0 :                CPABORT("gfn2-xtb tbd")
    1702              :             CASE DEFAULT
    1703         1228 :                CPABORT("GFN type")
    1704              :             END SELECT
    1705              :          END IF
    1706              :          CALL section_vals_val_get(xtb_parameter, "DISPERSION_PARAMETER_FILE", &
    1707         1228 :                                    c_val=qs_control%xtb_control%dispersion_parameter_file)
    1708              :          ! global parameters
    1709         1228 :          CALL section_vals_val_get(xtb_parameter, "HUCKEL_CONSTANTS", explicit=explicit)
    1710         1228 :          IF (explicit) THEN
    1711            0 :             CALL section_vals_val_get(xtb_parameter, "HUCKEL_CONSTANTS", r_vals=scal)
    1712            0 :             qs_control%xtb_control%ks = scal(1)
    1713            0 :             qs_control%xtb_control%kp = scal(2)
    1714            0 :             qs_control%xtb_control%kd = scal(3)
    1715            0 :             qs_control%xtb_control%ksp = scal(4)
    1716            0 :             qs_control%xtb_control%k2sh = scal(5)
    1717            0 :             IF (qs_control%xtb_control%gfn_type == 0) THEN
    1718              :                ! enforce ksp for gfn0
    1719            0 :                qs_control%xtb_control%ksp = 0.5_dp*(scal(1) + scal(2))
    1720              :             END IF
    1721              :          ELSE
    1722         1922 :             SELECT CASE (qs_control%xtb_control%gfn_type)
    1723              :             CASE (0)
    1724          694 :                qs_control%xtb_control%ks = 2.00_dp
    1725          694 :                qs_control%xtb_control%kp = 2.4868_dp
    1726          694 :                qs_control%xtb_control%kd = 2.27_dp
    1727          694 :                qs_control%xtb_control%ksp = 2.2434_dp
    1728          694 :                qs_control%xtb_control%k2sh = 1.1241_dp
    1729              :             CASE (1)
    1730          534 :                qs_control%xtb_control%ks = 1.85_dp
    1731          534 :                qs_control%xtb_control%kp = 2.25_dp
    1732          534 :                qs_control%xtb_control%kd = 2.00_dp
    1733          534 :                qs_control%xtb_control%ksp = 2.08_dp
    1734          534 :                qs_control%xtb_control%k2sh = 2.85_dp
    1735              :             CASE (2)
    1736            0 :                CPABORT("gfn2-xtb tbd")
    1737              :             CASE DEFAULT
    1738         1228 :                CPABORT("GFN type")
    1739              :             END SELECT
    1740              :          END IF
    1741         1228 :          CALL section_vals_val_get(xtb_parameter, "COULOMB_CONSTANTS", explicit=explicit)
    1742         1228 :          IF (explicit) THEN
    1743            0 :             CALL section_vals_val_get(xtb_parameter, "COULOMB_CONSTANTS", r_vals=scal)
    1744            0 :             qs_control%xtb_control%kg = scal(1)
    1745            0 :             qs_control%xtb_control%kf = scal(2)
    1746              :          ELSE
    1747         1922 :             SELECT CASE (qs_control%xtb_control%gfn_type)
    1748              :             CASE (0)
    1749          694 :                qs_control%xtb_control%kg = 2.00_dp
    1750          694 :                qs_control%xtb_control%kf = 1.50_dp
    1751              :             CASE (1)
    1752          534 :                qs_control%xtb_control%kg = 2.00_dp
    1753          534 :                qs_control%xtb_control%kf = 1.50_dp
    1754              :             CASE (2)
    1755            0 :                CPABORT("gfn2-xtb tbd")
    1756              :             CASE DEFAULT
    1757         1228 :                CPABORT("GFN type")
    1758              :             END SELECT
    1759              :          END IF
    1760         1228 :          CALL section_vals_val_get(xtb_parameter, "CN_CONSTANTS", r_vals=scal)
    1761         1228 :          qs_control%xtb_control%kcns = scal(1)
    1762         1228 :          qs_control%xtb_control%kcnp = scal(2)
    1763         1228 :          qs_control%xtb_control%kcnd = scal(3)
    1764              :          !
    1765         1228 :          CALL section_vals_val_get(xtb_parameter, "EN_CONSTANTS", explicit=explicit)
    1766         1228 :          IF (explicit) THEN
    1767            0 :             CALL section_vals_val_get(xtb_parameter, "EN_CONSTANTS", r_vals=scal)
    1768            0 :             SELECT CASE (qs_control%xtb_control%gfn_type)
    1769              :             CASE (0)
    1770            0 :                qs_control%xtb_control%ksen = scal(1)
    1771            0 :                qs_control%xtb_control%kpen = scal(2)
    1772            0 :                qs_control%xtb_control%kden = scal(3)
    1773              :             CASE (1)
    1774            0 :                qs_control%xtb_control%ken = scal(1)
    1775              :             CASE (2)
    1776            0 :                CPABORT("gfn2-xtb tbd")
    1777              :             CASE DEFAULT
    1778            0 :                CPABORT("GFN type")
    1779              :             END SELECT
    1780              :          ELSE
    1781         1922 :             SELECT CASE (qs_control%xtb_control%gfn_type)
    1782              :             CASE (0)
    1783          694 :                qs_control%xtb_control%ksen = 0.006_dp
    1784          694 :                qs_control%xtb_control%kpen = -0.001_dp
    1785          694 :                qs_control%xtb_control%kden = -0.002_dp
    1786              :             CASE (1)
    1787          534 :                qs_control%xtb_control%ken = -0.007_dp
    1788              :             CASE (2)
    1789            0 :                CPABORT("gfn2-xtb tbd")
    1790              :             CASE DEFAULT
    1791         1228 :                CPABORT("GFN type")
    1792              :             END SELECT
    1793              :          END IF
    1794              :          ! ben
    1795         1228 :          CALL section_vals_val_get(xtb_parameter, "BEN_CONSTANT", r_vals=scal)
    1796         1228 :          qs_control%xtb_control%ben = scal(1)
    1797              :          ! enscale (hidden parameter in repulsion
    1798         1228 :          CALL section_vals_val_get(xtb_parameter, "ENSCALE", explicit=explicit)
    1799         1228 :          IF (explicit) THEN
    1800              :             CALL section_vals_val_get(xtb_parameter, "ENSCALE", &
    1801            0 :                                       r_val=qs_control%xtb_control%enscale)
    1802              :          ELSE
    1803         1922 :             SELECT CASE (qs_control%xtb_control%gfn_type)
    1804              :             CASE (0)
    1805          694 :                qs_control%xtb_control%enscale = -0.09_dp
    1806              :             CASE (1)
    1807          534 :                qs_control%xtb_control%enscale = 0._dp
    1808              :             CASE (2)
    1809            0 :                CPABORT("gfn2-xtb tbd")
    1810              :             CASE DEFAULT
    1811         1228 :                CPABORT("GFN type")
    1812              :             END SELECT
    1813              :          END IF
    1814              :          ! XB
    1815              :          CALL section_vals_val_get(xtb_section, "USE_HALOGEN_CORRECTION", &
    1816         1228 :                                    l_val=qs_control%xtb_control%xb_interaction)
    1817         1228 :          CALL section_vals_val_get(xtb_parameter, "HALOGEN_BINDING", r_vals=scal)
    1818         1228 :          qs_control%xtb_control%kxr = scal(1)
    1819         1228 :          qs_control%xtb_control%kx2 = scal(2)
    1820              :          ! NONBONDED interactions
    1821              :          CALL section_vals_val_get(xtb_section, "DO_NONBONDED", &
    1822         1228 :                                    l_val=qs_control%xtb_control%do_nonbonded)
    1823         1228 :          CALL section_vals_get(nonbonded_section, explicit=explicit)
    1824         1228 :          IF (explicit .AND. qs_control%xtb_control%do_nonbonded) THEN
    1825            6 :             CALL section_vals_get(genpot_section, explicit=explicit, n_repetition=ngp)
    1826            6 :             IF (explicit) THEN
    1827            6 :                CALL pair_potential_reallocate(qs_control%xtb_control%nonbonded, 1, ngp, gp=.TRUE.)
    1828            6 :                CALL read_gp_section(qs_control%xtb_control%nonbonded, genpot_section, 0)
    1829              :             END IF
    1830              :          END IF !nonbonded
    1831              :          CALL section_vals_val_get(xtb_section, "EPS_PAIRPOTENTIAL", &
    1832         1228 :                                    r_val=qs_control%xtb_control%eps_pair)
    1833              :          ! SR Coulomb
    1834         1228 :          CALL section_vals_val_get(xtb_parameter, "COULOMB_SR_CUT", r_vals=scal)
    1835         1228 :          qs_control%xtb_control%coulomb_sr_cut = scal(1)
    1836         1228 :          CALL section_vals_val_get(xtb_parameter, "COULOMB_SR_EPS", r_vals=scal)
    1837         1228 :          qs_control%xtb_control%coulomb_sr_eps = scal(1)
    1838              :          ! XB_radius
    1839         1228 :          CALL section_vals_val_get(xtb_parameter, "XB_RADIUS", r_val=qs_control%xtb_control%xb_radius)
    1840              :          ! Kab
    1841         1228 :          CALL section_vals_val_get(xtb_parameter, "KAB_PARAM", n_rep_val=n_rep)
    1842              :          ! Coulomb
    1843         1922 :          SELECT CASE (qs_control%xtb_control%gfn_type)
    1844              :          CASE (0)
    1845          694 :             qs_control%xtb_control%coulomb_interaction = .FALSE.
    1846          694 :             qs_control%xtb_control%coulomb_lr = .FALSE.
    1847          694 :             qs_control%xtb_control%tb3_interaction = .FALSE.
    1848          694 :             qs_control%xtb_control%check_atomic_charges = .FALSE.
    1849              :             CALL section_vals_val_get(xtb_section, "VARIATIONAL_DIPOLE", &
    1850          694 :                                       l_val=qs_control%xtb_control%var_dipole)
    1851              :          CASE (1)
    1852              :             ! For debugging purposes
    1853              :             CALL section_vals_val_get(xtb_section, "COULOMB_INTERACTION", &
    1854          534 :                                       l_val=qs_control%xtb_control%coulomb_interaction)
    1855              :             CALL section_vals_val_get(xtb_section, "COULOMB_LR", &
    1856          534 :                                       l_val=qs_control%xtb_control%coulomb_lr)
    1857              :             CALL section_vals_val_get(xtb_section, "TB3_INTERACTION", &
    1858          534 :                                       l_val=qs_control%xtb_control%tb3_interaction)
    1859              :             ! Check for bad atomic charges
    1860              :             CALL section_vals_val_get(xtb_section, "CHECK_ATOMIC_CHARGES", &
    1861          534 :                                       l_val=qs_control%xtb_control%check_atomic_charges)
    1862          534 :             qs_control%xtb_control%var_dipole = .FALSE.
    1863              :          CASE (2)
    1864            0 :             CPABORT("gfn2-xtb tbd")
    1865              :          CASE DEFAULT
    1866         1228 :             CPABORT("GFN type")
    1867              :          END SELECT
    1868         1228 :          qs_control%xtb_control%kab_nval = n_rep
    1869         1228 :          IF (n_rep > 0) THEN
    1870            6 :             ALLOCATE (qs_control%xtb_control%kab_param(3, n_rep))
    1871            6 :             ALLOCATE (qs_control%xtb_control%kab_types(2, n_rep))
    1872            6 :             ALLOCATE (qs_control%xtb_control%kab_vals(n_rep))
    1873            4 :             DO j = 1, n_rep
    1874            2 :                CALL section_vals_val_get(xtb_parameter, "KAB_PARAM", i_rep_val=j, c_vals=clist)
    1875            2 :                qs_control%xtb_control%kab_param(1, j) = clist(1)
    1876              :                CALL get_ptable_info(clist(1) (1:2), &
    1877            2 :                                     ielement=qs_control%xtb_control%kab_types(1, j))
    1878            2 :                qs_control%xtb_control%kab_param(2, j) = clist(2)
    1879              :                CALL get_ptable_info(clist(2) (1:2), &
    1880            2 :                                     ielement=qs_control%xtb_control%kab_types(2, j))
    1881            2 :                qs_control%xtb_control%kab_param(3, j) = clist(3)
    1882            4 :                READ (clist(3), '(F10.0)') qs_control%xtb_control%kab_vals(j)
    1883              :             END DO
    1884              :          END IF
    1885              : 
    1886              :          ! Spin Polarisation
    1887         1228 :          CALL section_vals_val_get(xtb_parameter, "SPIN_POL_PARAM", n_rep_val=n_rep)
    1888         1228 :          IF (n_rep > 0) THEN
    1889            6 :             ALLOCATE (qs_control%xtb_control%spinpol_type(n_rep))
    1890            6 :             ALLOCATE (qs_control%xtb_control%spinpol_vals(6, n_rep))
    1891            8 :             DO j = 1, n_rep
    1892            6 :                CALL section_vals_val_get(xtb_parameter, "SPIN_POL_PARAM", i_rep_val=j, c_vals=clist)
    1893            6 :                READ (clist(1), '(A)') cval
    1894            6 :                element_symbol = ADJUSTL(TRIM(cval))
    1895            6 :                CALL get_ptable_info(element_symbol, znum)
    1896            6 :                qs_control%xtb_control%spinpol_type(j) = znum
    1897            6 :                READ (clist(2), '(F20.8)') qs_control%xtb_control%spinpol_vals(1, j)
    1898            6 :                READ (clist(3), '(F20.8)') qs_control%xtb_control%spinpol_vals(2, j)
    1899            6 :                READ (clist(4), '(F20.8)') qs_control%xtb_control%spinpol_vals(3, j)
    1900            6 :                READ (clist(5), '(F20.8)') qs_control%xtb_control%spinpol_vals(4, j)
    1901            6 :                READ (clist(6), '(F20.8)') qs_control%xtb_control%spinpol_vals(5, j)
    1902           14 :                READ (clist(7), '(F20.8)') qs_control%xtb_control%spinpol_vals(6, j)
    1903              :             END DO
    1904              :          END IF
    1905              : 
    1906         1228 :          IF (qs_control%xtb_control%gfn_type == 0) THEN
    1907          694 :             CALL section_vals_val_get(xtb_parameter, "SRB_PARAMETER", r_vals=scal)
    1908          694 :             qs_control%xtb_control%ksrb = scal(1)
    1909          694 :             qs_control%xtb_control%esrb = scal(2)
    1910          694 :             qs_control%xtb_control%gscal = scal(3)
    1911          694 :             qs_control%xtb_control%c1srb = scal(4)
    1912          694 :             qs_control%xtb_control%c2srb = scal(5)
    1913          694 :             qs_control%xtb_control%shift = scal(6)
    1914              :          END IF
    1915              : 
    1916         1228 :          CALL section_vals_val_get(xtb_section, "EN_SHIFT_TYPE", c_val=cval)
    1917         1228 :          CALL uppercase(cval)
    1918         1228 :          SELECT CASE (TRIM(cval))
    1919              :          CASE ("SELECT")
    1920            0 :             qs_control%xtb_control%enshift_type = 0
    1921              :          CASE ("MOLECULE")
    1922         1228 :             qs_control%xtb_control%enshift_type = 1
    1923              :          CASE ("CRYSTAL")
    1924            0 :             qs_control%xtb_control%enshift_type = 2
    1925              :          CASE DEFAULT
    1926         1228 :             CPABORT("Unknown value for EN_SHIFT_TYPE")
    1927              :          END SELECT
    1928              : 
    1929              :          ! EEQ solver params
    1930         1228 :          CALL read_eeq_param(eeq_section, qs_control%xtb_control%eeq_sparam)
    1931              :       END IF
    1932              : 
    1933              :       ! Optimize LRI basis set
    1934         9032 :       CALL section_vals_get(lri_optbas_section, explicit=qs_control%lri_optbas)
    1935              : 
    1936              :       ! Use tblite if selected through XTB/GFN_TYPE TBLITE.
    1937         9032 :       IF (qs_control%xtb_control%do_tblite) THEN
    1938              :          CALL section_vals_val_get(xtb_tblite, "METHOD", &
    1939          188 :                                    i_val=qs_control%xtb_control%tblite_method)
    1940              :          CALL section_vals_val_get(xtb_tblite, "PARAM", &
    1941          188 :                                    c_val=qs_control%xtb_control%tblite_param_file)
    1942              :          CALL section_vals_val_get(xtb_tblite, "ACCURACY", &
    1943          188 :                                    r_val=qs_control%xtb_control%tblite_accuracy)
    1944          188 :          IF (qs_control%xtb_control%tblite_accuracy <= 0.0_dp) THEN
    1945            0 :             CPABORT("XTB/TBLITE/ACCURACY must be positive")
    1946              :          END IF
    1947          188 :          IF (qs_control%xtb_control%tblite_mixer_damping <= 0.0_dp) THEN
    1948            0 :             CPABORT("XTB/TBLITE_MIXER/DAMPING must be positive")
    1949              :          END IF
    1950          188 :          CALL section_vals_val_get(xtb_tblite, "REFERENCE_CLI", l_val=tblite_reference_cli)
    1951          188 :          CALL section_vals_get(xtb_tblite_ref_cli, explicit=tblite_reference_cli_section)
    1952          188 :          IF (tblite_reference_cli .AND. (.NOT. tblite_reference_cli_section)) THEN
    1953            0 :             CPABORT("XTB/TBLITE/REFERENCE_CLI keyword requires an XTB/TBLITE/REFERENCE_CLI section")
    1954              :          END IF
    1955          188 :          IF (tblite_reference_cli .OR. tblite_reference_cli_section) THEN
    1956            2 :             CALL read_xtb_reference_cli_section(xtb_tblite_ref_cli, qs_control%xtb_control%reference_cli, cell)
    1957            2 :             qs_control%xtb_control%reference_cli%enabled = .TRUE.
    1958              :          END IF
    1959          188 :          CALL cite_reference(Katbashev2025)
    1960              :          ! tblite handles periodic long-range terms internally from the CP2K cell periodicity.
    1961              :          ! Keep xtb_control%do_ewald as read above from XTB/DO_EWALD or SUBSYS/CELL/PERIODIC,
    1962              :          ! matching the DFTB and CP2K-internal xTB setup.
    1963              :       END IF
    1964              : 
    1965         9032 :       CALL timestop(handle)
    1966         9032 :    END SUBROUTINE read_qs_section
    1967              : 
    1968              : ! **************************************************************************************************
    1969              : !> \brief Read a TBLITE_MIXER section.
    1970              : !> \param mixer_section input section
    1971              : !> \param iterations tblite SCC iteration limit
    1972              : !> \param memory Broyden history length
    1973              : !> \param solver native tblite electronic solver id
    1974              : !> \param damping mixer damping parameter
    1975              : !> \param omega0 Broyden regularization weight
    1976              : !> \param min_weight minimum dynamic Broyden weight
    1977              : !> \param max_weight maximum dynamic Broyden weight
    1978              : !> \param weight_factor residual-to-weight scaling factor
    1979              : !> \param section_name diagnostic section name
    1980              : ! **************************************************************************************************
    1981         1558 :    SUBROUTINE read_tblite_mixer_section(mixer_section, iterations, memory, solver, damping, omega0, min_weight, &
    1982              :                                         max_weight, weight_factor, section_name)
    1983              :       TYPE(section_vals_type), POINTER                   :: mixer_section
    1984              :       INTEGER, INTENT(INOUT)                             :: iterations, memory, solver
    1985              :       REAL(KIND=dp), INTENT(INOUT)                       :: damping, omega0, min_weight, max_weight, &
    1986              :                                                             weight_factor
    1987              :       CHARACTER(LEN=*), INTENT(IN)                       :: section_name
    1988              : 
    1989              :       LOGICAL                                            :: explicit, memory_explicit
    1990              : 
    1991         1526 :       CALL section_vals_get(mixer_section, explicit=explicit)
    1992         1526 :       IF (.NOT. explicit) RETURN
    1993              : 
    1994            4 :       CALL section_vals_val_get(mixer_section, "ITERATIONS", i_val=iterations)
    1995            4 :       CALL section_vals_val_get(mixer_section, "MEMORY", explicit=memory_explicit)
    1996            4 :       IF (memory_explicit) CALL section_vals_val_get(mixer_section, "MEMORY", i_val=memory)
    1997            4 :       IF (.NOT. memory_explicit .OR. memory == tblite_mixer_memory_inherit) THEN
    1998            2 :          memory = iterations
    1999              :       END IF
    2000            4 :       CALL section_vals_val_get(mixer_section, "SOLVER", i_val=solver)
    2001            4 :       CALL section_vals_val_get(mixer_section, "DAMPING", r_val=damping)
    2002            4 :       CALL section_vals_val_get(mixer_section, "OMEGA0", r_val=omega0)
    2003            4 :       CALL section_vals_val_get(mixer_section, "MIN_WEIGHT", r_val=min_weight)
    2004            4 :       CALL section_vals_val_get(mixer_section, "MAX_WEIGHT", r_val=max_weight)
    2005            4 :       CALL section_vals_val_get(mixer_section, "WEIGHT_FACTOR", r_val=weight_factor)
    2006              : 
    2007            4 :       IF (iterations < 1) CPABORT(TRIM(section_name)//"/ITERATIONS must be positive")
    2008            4 :       IF (memory < 1) CPABORT(TRIM(section_name)//"/MEMORY must be positive")
    2009            4 :       SELECT CASE (solver)
    2010              :       CASE (tblite_solver_gvd, tblite_solver_gvr)
    2011              :       CASE DEFAULT
    2012            4 :          CPABORT(TRIM(section_name)//"/SOLVER must be GVD or GVR")
    2013              :       END SELECT
    2014            4 :       IF (damping <= 0.0_dp) CPABORT(TRIM(section_name)//"/DAMPING must be positive")
    2015            4 :       IF (omega0 <= 0.0_dp) CPABORT(TRIM(section_name)//"/OMEGA0 must be positive")
    2016            4 :       IF (min_weight <= 0.0_dp) CPABORT(TRIM(section_name)//"/MIN_WEIGHT must be positive")
    2017            4 :       IF (max_weight <= 0.0_dp) CPABORT(TRIM(section_name)//"/MAX_WEIGHT must be positive")
    2018            4 :       IF (max_weight < min_weight) THEN
    2019            0 :          CPABORT(TRIM(section_name)//"/MAX_WEIGHT must not be smaller than MIN_WEIGHT")
    2020              :       END IF
    2021            4 :       IF (weight_factor <= 0.0_dp) CPABORT(TRIM(section_name)//"/WEIGHT_FACTOR must be positive")
    2022              : 
    2023              :    END SUBROUTINE read_tblite_mixer_section
    2024              : 
    2025              : ! **************************************************************************************************
    2026              : !> \brief Read native tblite CLI reference options.
    2027              : !> \param ref_cli_section input section
    2028              : !> \param ref_cli reference CLI control data
    2029              : !> \param cell optional cell used to transform Cartesian input vectors
    2030              : ! **************************************************************************************************
    2031            2 :    SUBROUTINE read_xtb_reference_cli_section(ref_cli_section, ref_cli, cell)
    2032              :       TYPE(section_vals_type), POINTER                   :: ref_cli_section
    2033              :       TYPE(xtb_reference_cli_type), INTENT(INOUT)        :: ref_cli
    2034              :       TYPE(cell_type), OPTIONAL, POINTER                 :: cell
    2035              : 
    2036            2 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: efield
    2037              :       TYPE(section_vals_type), POINTER                   :: fit_section, guess_section, &
    2038              :                                                             param_section, solvation_section, &
    2039              :                                                             tagdiff_section
    2040              : 
    2041            2 :       CALL section_vals_val_get(ref_cli_section, "_SECTION_PARAMETERS_", l_val=ref_cli%enabled)
    2042            2 :       CALL section_vals_val_get(ref_cli_section, "PROGRAM_NAME", c_val=ref_cli%program_name)
    2043            2 :       CALL section_vals_val_get(ref_cli_section, "GUESS", i_val=ref_cli%guess)
    2044            2 :       CALL section_vals_val_get(ref_cli_section, "WORK_DIRECTORY", c_val=ref_cli%work_directory)
    2045            2 :       CALL section_vals_val_get(ref_cli_section, "PREFIX", c_val=ref_cli%prefix)
    2046            2 :       CALL section_vals_val_get(ref_cli_section, "INPUT_FORMAT", c_val=ref_cli%input_format)
    2047            2 :       CALL section_vals_val_get(ref_cli_section, "RESTART", c_val=ref_cli%restart_file)
    2048            2 :       CALL section_vals_val_get(ref_cli_section, "GRAD", c_val=ref_cli%grad_file)
    2049            2 :       CALL section_vals_val_get(ref_cli_section, "JSON", c_val=ref_cli%json_file)
    2050            2 :       CALL section_vals_val_get(ref_cli_section, "POST_PROCESSING", c_val=ref_cli%post_processing)
    2051              :       CALL section_vals_val_get(ref_cli_section, "POST_PROCESSING_OUTPUT", &
    2052            2 :                                 c_val=ref_cli%post_processing_output_file)
    2053            2 :       CALL section_vals_val_get(ref_cli_section, "EFIELD", explicit=ref_cli%efield_active)
    2054            2 :       IF (ref_cli%efield_active) THEN
    2055            2 :          NULLIFY (efield)
    2056            2 :          CALL section_vals_val_get(ref_cli_section, "EFIELD", r_vals=efield)
    2057            8 :          ref_cli%efield = efield(1:3)
    2058            2 :          IF (PRESENT(cell)) THEN
    2059            2 :             IF (ASSOCIATED(cell)) CALL cell_transform_input_cartesian(cell, ref_cli%efield)
    2060              :          END IF
    2061              :       END IF
    2062            2 :       solvation_section => section_vals_get_subs_vals(ref_cli_section, "IMPLICIT_SOLVATION")
    2063            2 :       CALL section_vals_get(solvation_section, explicit=ref_cli%solvation_active)
    2064            2 :       IF (ref_cli%solvation_active) THEN
    2065            2 :          CALL section_vals_val_get(solvation_section, "MODEL", i_val=ref_cli%solvation_model)
    2066            2 :          CALL section_vals_val_get(solvation_section, "SOLVENT", c_val=ref_cli%solvation_solvent)
    2067            2 :          CALL section_vals_val_get(solvation_section, "BORN_KERNEL", i_val=ref_cli%solvation_born_kernel)
    2068            2 :          CALL section_vals_val_get(solvation_section, "SOLUTION_STATE", i_val=ref_cli%solvation_state)
    2069            2 :          IF (LEN_TRIM(ref_cli%solvation_solvent) == 0) THEN
    2070            0 :             CPABORT("REFERENCE_CLI implicit solvation needs SOLVENT")
    2071              :          END IF
    2072            2 :          IF (ref_cli%solvation_model == tblite_cli_solvation_cpcm .AND. &
    2073              :              ref_cli%solvation_born_kernel /= tblite_cli_born_kernel_auto) THEN
    2074            0 :             CPABORT("BORN_KERNEL is invalid with MODEL CPCM")
    2075              :          END IF
    2076            2 :          IF (ref_cli%solvation_state /= tblite_cli_solution_state_gsolv) THEN
    2077            2 :             SELECT CASE (ref_cli%solvation_model)
    2078              :             CASE (tblite_cli_solvation_alpb, tblite_cli_solvation_gbsa)
    2079              :                ! Native tblite supports solution-state shifts for parametrized named-solvent ALPB/GBSA.
    2080              :             CASE (tblite_cli_solvation_gbe, tblite_cli_solvation_gb, tblite_cli_solvation_cpcm)
    2081            2 :                CPABORT("SOLUTION_STATE is valid only for ALPB/GBSA")
    2082              :             END SELECT
    2083              :          END IF
    2084              :       END IF
    2085              :       CALL section_vals_val_get(ref_cli_section, "ELECTRONIC_TEMPERATURE_GUESS", &
    2086            2 :                                 r_val=ref_cli%electronic_temperature_guess)
    2087            2 :       IF (ref_cli%electronic_temperature_guess < 0.0_dp) THEN
    2088            0 :          CPABORT("XTB/TBLITE/REFERENCE_CLI/ELECTRONIC_TEMPERATURE_GUESS must not be negative")
    2089              :       END IF
    2090            2 :       IF (ref_cli%electronic_temperature_guess > 0.0_dp .AND. ref_cli%guess /= tblite_guess_ceh) THEN
    2091            0 :          CPABORT("XTB/TBLITE/REFERENCE_CLI/ELECTRONIC_TEMPERATURE_GUESS requires GUESS CEH")
    2092              :       END IF
    2093            2 :       guess_section => section_vals_get_subs_vals(ref_cli_section, "GUESS_CLI")
    2094            2 :       CALL section_vals_get(guess_section, explicit=ref_cli%guess_cli%enabled)
    2095            2 :       IF (ref_cli%guess_cli%enabled) THEN
    2096            0 :          CALL section_vals_val_get(guess_section, "METHOD", i_val=ref_cli%guess_cli%method)
    2097              :          CALL section_vals_val_get(guess_section, "ELECTRONIC_TEMPERATURE_GUESS", &
    2098            0 :                                    r_val=ref_cli%guess_cli%electronic_temperature_guess)
    2099            0 :          IF (ref_cli%guess_cli%electronic_temperature_guess < 0.0_dp) THEN
    2100            0 :             CPABORT("REFERENCE_CLI/GUESS_CLI/ELECTRONIC_TEMPERATURE_GUESS must not be negative")
    2101              :          END IF
    2102            0 :          CALL section_vals_val_get(guess_section, "SOLVER", i_val=ref_cli%guess_cli%solver)
    2103            0 :          CALL section_vals_val_get(guess_section, "EFIELD", explicit=ref_cli%guess_cli%efield_active)
    2104            0 :          IF (ref_cli%guess_cli%efield_active) THEN
    2105            0 :             NULLIFY (efield)
    2106            0 :             CALL section_vals_val_get(guess_section, "EFIELD", r_vals=efield)
    2107            0 :             ref_cli%guess_cli%efield = efield(1:3)
    2108            0 :             IF (PRESENT(cell)) THEN
    2109            0 :                IF (ASSOCIATED(cell)) CALL cell_transform_input_cartesian(cell, ref_cli%guess_cli%efield)
    2110              :             END IF
    2111              :          END IF
    2112            0 :          CALL section_vals_val_get(guess_section, "GRAD", l_val=ref_cli%guess_cli%grad)
    2113            0 :          CALL section_vals_val_get(guess_section, "JSON", c_val=ref_cli%guess_cli%json_file)
    2114            0 :          CALL section_vals_val_get(guess_section, "INPUT_FORMAT", c_val=ref_cli%guess_cli%input_format)
    2115            0 :          CALL section_vals_val_get(guess_section, "INPUT_FILE", c_val=ref_cli%guess_cli%input_file)
    2116              :       END IF
    2117            2 :       param_section => section_vals_get_subs_vals(ref_cli_section, "PARAM_CLI")
    2118            2 :       CALL section_vals_get(param_section, explicit=ref_cli%param_cli%enabled)
    2119            2 :       IF (ref_cli%param_cli%enabled) THEN
    2120              :          CALL section_vals_val_get(param_section, "METHOD", explicit=ref_cli%param_cli%method_explicit, &
    2121            0 :                                    i_val=ref_cli%param_cli%method)
    2122            0 :          CALL section_vals_val_get(param_section, "OUTPUT", c_val=ref_cli%param_cli%output_file)
    2123            0 :          CALL section_vals_val_get(param_section, "INPUT_FILE", c_val=ref_cli%param_cli%input_file)
    2124              :       END IF
    2125            2 :       fit_section => section_vals_get_subs_vals(ref_cli_section, "FIT_CLI")
    2126            2 :       CALL section_vals_get(fit_section, explicit=ref_cli%fit_cli%enabled)
    2127            2 :       IF (ref_cli%fit_cli%enabled) THEN
    2128            0 :          CALL section_vals_val_get(fit_section, "PARAM_FILE", c_val=ref_cli%fit_cli%param_file)
    2129            0 :          CALL section_vals_val_get(fit_section, "INPUT_FILE", c_val=ref_cli%fit_cli%input_file)
    2130            0 :          CALL section_vals_val_get(fit_section, "DRY_RUN", l_val=ref_cli%fit_cli%dry_run)
    2131            0 :          CALL section_vals_val_get(fit_section, "COPY", c_val=ref_cli%fit_cli%copy_file)
    2132            0 :          IF (LEN_TRIM(ref_cli%fit_cli%param_file) == 0) THEN
    2133            0 :             CPABORT("XTB/TBLITE/REFERENCE_CLI/FIT_CLI needs PARAM_FILE")
    2134              :          END IF
    2135            0 :          IF (LEN_TRIM(ref_cli%fit_cli%input_file) == 0) THEN
    2136            0 :             CPABORT("XTB/TBLITE/REFERENCE_CLI/FIT_CLI needs INPUT_FILE")
    2137              :          END IF
    2138              :       END IF
    2139            2 :       tagdiff_section => section_vals_get_subs_vals(ref_cli_section, "TAGDIFF_CLI")
    2140            2 :       CALL section_vals_get(tagdiff_section, explicit=ref_cli%tagdiff_cli%enabled)
    2141            2 :       IF (ref_cli%tagdiff_cli%enabled) THEN
    2142            0 :          CALL section_vals_val_get(tagdiff_section, "ACTUAL", c_val=ref_cli%tagdiff_cli%actual_file)
    2143            0 :          CALL section_vals_val_get(tagdiff_section, "REFERENCE", c_val=ref_cli%tagdiff_cli%reference_file)
    2144            0 :          CALL section_vals_val_get(tagdiff_section, "FIT", l_val=ref_cli%tagdiff_cli%fit)
    2145            0 :          IF (LEN_TRIM(ref_cli%tagdiff_cli%actual_file) == 0) THEN
    2146            0 :             CPABORT("XTB/TBLITE/REFERENCE_CLI/TAGDIFF_CLI needs ACTUAL")
    2147              :          END IF
    2148            0 :          IF (LEN_TRIM(ref_cli%tagdiff_cli%reference_file) == 0) THEN
    2149            0 :             CPABORT("XTB/TBLITE/REFERENCE_CLI/TAGDIFF_CLI needs REFERENCE")
    2150              :          END IF
    2151              :       END IF
    2152            2 :       CALL section_vals_val_get(ref_cli_section, "KEEP_FILES", l_val=ref_cli%keep_files)
    2153            2 :       CALL section_vals_val_get(ref_cli_section, "ERROR_LIMIT", r_val=ref_cli%error_limit)
    2154            2 :       CALL section_vals_val_get(ref_cli_section, "STOP_ON_ERROR", l_val=ref_cli%stop_on_error)
    2155            2 :       CALL section_vals_val_get(ref_cli_section, "CHECK_ENERGY", l_val=ref_cli%check_energy)
    2156            2 :       CALL section_vals_val_get(ref_cli_section, "CHECK_FORCES", l_val=ref_cli%check_forces)
    2157            2 :       CALL section_vals_val_get(ref_cli_section, "CHECK_VIRIAL", l_val=ref_cli%check_virial)
    2158              : 
    2159            2 :    END SUBROUTINE read_xtb_reference_cli_section
    2160              : 
    2161              : ! **************************************************************************************************
    2162              : !> \brief Read TDDFPT-related input parameters.
    2163              : !> \param t_control  TDDFPT control parameters
    2164              : !> \param t_section  TDDFPT input section
    2165              : !> \param qs_control Quickstep control parameters
    2166              : ! **************************************************************************************************
    2167         9064 :    SUBROUTINE read_tddfpt2_control(t_control, t_section, qs_control)
    2168              :       TYPE(tddfpt2_control_type), POINTER                :: t_control
    2169              :       TYPE(section_vals_type), POINTER                   :: t_section
    2170              :       TYPE(qs_control_type), POINTER                     :: qs_control
    2171              : 
    2172              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'read_tddfpt2_control'
    2173              : 
    2174              :       CHARACTER(LEN=default_string_length), &
    2175         9064 :          DIMENSION(:), POINTER                           :: tmpstringlist
    2176              :       INTEGER                                            :: handle, irep, isize, nrep
    2177         9064 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: inds
    2178              :       LOGICAL                                            :: do_ewald, do_exchange, expl, explicit, &
    2179              :                                                             multigrid_set
    2180              :       REAL(KIND=dp)                                      :: filter, fval, hfx
    2181              :       TYPE(section_vals_type), POINTER                   :: dipole_section, mgrid_section, &
    2182              :                                                             soc_section, stda_section, xc_func, &
    2183              :                                                             xc_section
    2184              : 
    2185         9064 :       CALL timeset(routineN, handle)
    2186              : 
    2187         9064 :       CALL section_vals_val_get(t_section, "_SECTION_PARAMETERS_", l_val=t_control%enabled)
    2188              : 
    2189         9064 :       CALL section_vals_val_get(t_section, "NSTATES", i_val=t_control%nstates)
    2190         9064 :       CALL section_vals_val_get(t_section, "MAX_ITER", i_val=t_control%niters)
    2191         9064 :       CALL section_vals_val_get(t_section, "MAX_KV", i_val=t_control%nkvs)
    2192         9064 :       CALL section_vals_val_get(t_section, "NLUMO", i_val=t_control%nlumo)
    2193         9064 :       CALL section_vals_val_get(t_section, "NPROC_STATE", i_val=t_control%nprocs)
    2194         9064 :       CALL section_vals_val_get(t_section, "KERNEL", i_val=t_control%kernel)
    2195         9064 :       CALL section_vals_val_get(t_section, "SPINFLIP", i_val=t_control%spinflip)
    2196         9064 :       CALL section_vals_val_get(t_section, "OE_CORR", i_val=t_control%oe_corr)
    2197         9064 :       CALL section_vals_val_get(t_section, "EV_SHIFT", r_val=t_control%ev_shift)
    2198         9064 :       CALL section_vals_val_get(t_section, "EOS_SHIFT", r_val=t_control%eos_shift)
    2199              : 
    2200         9064 :       CALL section_vals_val_get(t_section, "CONVERGENCE", r_val=t_control%conv)
    2201         9064 :       CALL section_vals_val_get(t_section, "MIN_AMPLITUDE", r_val=t_control%min_excitation_amplitude)
    2202         9064 :       CALL section_vals_val_get(t_section, "ORTHOGONAL_EPS", r_val=t_control%orthogonal_eps)
    2203              : 
    2204         9064 :       CALL section_vals_val_get(t_section, "RESTART", l_val=t_control%is_restart)
    2205         9064 :       CALL section_vals_val_get(t_section, "RKS_TRIPLETS", l_val=t_control%rks_triplets)
    2206         9064 :       CALL section_vals_val_get(t_section, "DO_LRIGPW", l_val=t_control%do_lrigpw)
    2207         9064 :       CALL section_vals_val_get(t_section, "DO_SMEARING", l_val=t_control%do_smearing)
    2208         9064 :       CALL section_vals_val_get(t_section, "DO_BSE", l_val=t_control%do_bse)
    2209         9064 :       CALL section_vals_val_get(t_section, "DO_BSE_W_ONLY", l_val=t_control%do_bse_w_only)
    2210         9064 :       CALL section_vals_val_get(t_section, "DO_BSE_GW_ONLY", l_val=t_control%do_bse_gw_only)
    2211         9064 :       CALL section_vals_val_get(t_section, "ADMM_KERNEL_CORRECTION_SYMMETRIC", l_val=t_control%admm_symm)
    2212         9064 :       CALL section_vals_val_get(t_section, "ADMM_KERNEL_XC_CORRECTION", l_val=t_control%admm_xc_correction)
    2213         9064 :       CALL section_vals_val_get(t_section, "EXCITON_DESCRIPTORS", l_val=t_control%do_exciton_descriptors)
    2214         9064 :       CALL section_vals_val_get(t_section, "DIRECTIONAL_EXCITON_DESCRIPTORS", l_val=t_control%do_directional_exciton_descriptors)
    2215              : 
    2216              :       ! read automatically generated auxiliary basis for LRI
    2217         9064 :       CALL section_vals_val_get(t_section, "AUTO_BASIS", n_rep_val=nrep)
    2218        18128 :       DO irep = 1, nrep
    2219         9064 :          CALL section_vals_val_get(t_section, "AUTO_BASIS", i_rep_val=irep, c_vals=tmpstringlist)
    2220        18128 :          IF (SIZE(tmpstringlist) == 2) THEN
    2221         9064 :             CALL uppercase(tmpstringlist(2))
    2222        18128 :             SELECT CASE (tmpstringlist(2))
    2223              :             CASE ("X")
    2224         9064 :                SELECT CASE (tmpstringlist(1))
    2225              :                CASE ("X")
    2226              :                   ! Do nothing
    2227              :                CASE DEFAULT
    2228              :                   CALL cp_abort(__LOCATION__, &
    2229              :                                 "AUTO_BASIS: the size <X> is invalid for the "// &
    2230              :                                 "type <"//TRIM(ADJUSTL(tmpstringlist(1)))//">; "// &
    2231              :                                 "use one of SMALL, MEDIUM, LARGE, HUGE for "// &
    2232              :                                 "the size. The syntax AUTO_BASIS X X is a "// &
    2233              :                                 "reserved case for using NO automatically "// &
    2234         9064 :                                 "generated basis sets.")
    2235              :                END SELECT
    2236              :             CASE ("SMALL")
    2237            0 :                isize = 0
    2238              :             CASE ("MEDIUM")
    2239            0 :                isize = 1
    2240              :             CASE ("LARGE")
    2241            0 :                isize = 2
    2242              :             CASE ("HUGE")
    2243            0 :                isize = 3
    2244              :             CASE DEFAULT
    2245         9064 :                CPABORT("Unknown basis size in AUTO_BASIS keyword:"//TRIM(tmpstringlist(1)))
    2246              :             END SELECT
    2247              :             !
    2248         9064 :             SELECT CASE (tmpstringlist(1))
    2249              :             CASE ("X")
    2250              :             CASE ("P_LRI_AUX")
    2251            0 :                t_control%auto_basis_p_lri_aux = isize
    2252              :             CASE DEFAULT
    2253         9064 :                CPABORT("Unknown basis type in AUTO_BASIS keyword:"//TRIM(tmpstringlist(1)))
    2254              :             END SELECT
    2255              :          ELSE
    2256              :             CALL cp_abort(__LOCATION__, &
    2257            0 :                           "AUTO_BASIS keyword in &PROPERTIES &TDDFT section has a wrong number of arguments.")
    2258              :          END IF
    2259              :       END DO
    2260              : 
    2261         9064 :       IF (t_control%conv < 0) THEN
    2262            0 :          t_control%conv = ABS(t_control%conv)
    2263              :       END IF
    2264              : 
    2265              :       ! DIPOLE_MOMENTS subsection
    2266         9064 :       dipole_section => section_vals_get_subs_vals(t_section, "DIPOLE_MOMENTS")
    2267         9064 :       CALL section_vals_val_get(dipole_section, "DIPOLE_FORM", explicit=explicit)
    2268         9064 :       IF (explicit) THEN
    2269           36 :          CALL section_vals_val_get(dipole_section, "DIPOLE_FORM", i_val=t_control%dipole_form)
    2270              :       ELSE
    2271         9028 :          t_control%dipole_form = 0
    2272              :       END IF
    2273         9064 :       CALL section_vals_val_get(dipole_section, "REFERENCE", i_val=t_control%dipole_reference)
    2274         9064 :       CALL section_vals_val_get(dipole_section, "REFERENCE_POINT", explicit=explicit)
    2275         9064 :       IF (explicit) THEN
    2276            0 :          CALL section_vals_val_get(dipole_section, "REFERENCE_POINT", r_vals=t_control%dipole_ref_point)
    2277              :       ELSE
    2278         9064 :          NULLIFY (t_control%dipole_ref_point)
    2279         9064 :          IF (t_control%dipole_form == tddfpt_dipole_length .AND. t_control%dipole_reference == use_mom_ref_user) THEN
    2280            0 :             CPABORT("User-defined reference point should be given explicitly")
    2281              :          END IF
    2282              :       END IF
    2283              : 
    2284              :       !SOC subsection
    2285         9064 :       soc_section => section_vals_get_subs_vals(t_section, "SOC")
    2286         9064 :       CALL section_vals_get(soc_section, explicit=explicit)
    2287         9064 :       IF (explicit) THEN
    2288           10 :          t_control%do_soc = .TRUE.
    2289              :       END IF
    2290              : 
    2291              :       ! MGRID subsection
    2292         9064 :       mgrid_section => section_vals_get_subs_vals(t_section, "MGRID")
    2293         9064 :       CALL section_vals_get(mgrid_section, explicit=t_control%mgrid_is_explicit)
    2294              : 
    2295         9064 :       IF (t_control%mgrid_is_explicit) THEN
    2296           10 :          CALL section_vals_val_get(mgrid_section, "NGRIDS", i_val=t_control%mgrid_ngrids, explicit=explicit)
    2297           10 :          IF (.NOT. explicit) t_control%mgrid_ngrids = SIZE(qs_control%e_cutoff)
    2298              : 
    2299           10 :          CALL section_vals_val_get(mgrid_section, "CUTOFF", r_val=t_control%mgrid_cutoff, explicit=explicit)
    2300           10 :          IF (.NOT. explicit) t_control%mgrid_cutoff = qs_control%cutoff
    2301              : 
    2302              :          CALL section_vals_val_get(mgrid_section, "PROGRESSION_FACTOR", &
    2303           10 :                                    r_val=t_control%mgrid_progression_factor, explicit=explicit)
    2304           10 :          IF (explicit) THEN
    2305            0 :             IF (t_control%mgrid_progression_factor <= 1.0_dp) THEN
    2306              :                CALL cp_abort(__LOCATION__, &
    2307            0 :                              "Progression factor should be greater then 1.0 to ensure multi-grid ordering")
    2308              :             END IF
    2309              :          ELSE
    2310           10 :             t_control%mgrid_progression_factor = qs_control%progression_factor
    2311              :          END IF
    2312              : 
    2313           10 :          CALL section_vals_val_get(mgrid_section, "COMMENSURATE", l_val=t_control%mgrid_commensurate_mgrids, explicit=explicit)
    2314           10 :          IF (.NOT. explicit) t_control%mgrid_commensurate_mgrids = qs_control%commensurate_mgrids
    2315           10 :          IF (t_control%mgrid_commensurate_mgrids) THEN
    2316            0 :             IF (explicit) THEN
    2317            0 :                t_control%mgrid_progression_factor = 4.0_dp
    2318              :             ELSE
    2319            0 :                t_control%mgrid_progression_factor = qs_control%progression_factor
    2320              :             END IF
    2321              :          END IF
    2322              : 
    2323           10 :          CALL section_vals_val_get(mgrid_section, "REL_CUTOFF", r_val=t_control%mgrid_relative_cutoff, explicit=explicit)
    2324           10 :          IF (.NOT. explicit) t_control%mgrid_relative_cutoff = qs_control%relative_cutoff
    2325              : 
    2326           10 :          CALL section_vals_val_get(mgrid_section, "MULTIGRID_SET", l_val=multigrid_set, explicit=explicit)
    2327           10 :          IF (.NOT. explicit) multigrid_set = .FALSE.
    2328           10 :          IF (multigrid_set) THEN
    2329            0 :             CALL section_vals_val_get(mgrid_section, "MULTIGRID_CUTOFF", r_vals=t_control%mgrid_e_cutoff)
    2330              :          ELSE
    2331           10 :             NULLIFY (t_control%mgrid_e_cutoff)
    2332              :          END IF
    2333              : 
    2334           10 :          CALL section_vals_val_get(mgrid_section, "REALSPACE", l_val=t_control%mgrid_realspace_mgrids, explicit=explicit)
    2335           10 :          IF (.NOT. explicit) t_control%mgrid_realspace_mgrids = qs_control%realspace_mgrids
    2336              : 
    2337              :          CALL section_vals_val_get(mgrid_section, "SKIP_LOAD_BALANCE_DISTRIBUTED", &
    2338           10 :                                    l_val=t_control%mgrid_skip_load_balance, explicit=explicit)
    2339           10 :          IF (.NOT. explicit) t_control%mgrid_skip_load_balance = qs_control%skip_load_balance_distributed
    2340              : 
    2341           10 :          IF (ASSOCIATED(t_control%mgrid_e_cutoff)) THEN
    2342            0 :             IF (SIZE(t_control%mgrid_e_cutoff) /= t_control%mgrid_ngrids) THEN
    2343            0 :                CPABORT("Inconsistent values for number of multi-grids")
    2344              :             END IF
    2345              : 
    2346              :             ! sort multi-grids in descending order according to their cutoff values
    2347            0 :             t_control%mgrid_e_cutoff = -t_control%mgrid_e_cutoff
    2348            0 :             ALLOCATE (inds(t_control%mgrid_ngrids))
    2349            0 :             CALL sort(t_control%mgrid_e_cutoff, t_control%mgrid_ngrids, inds)
    2350            0 :             DEALLOCATE (inds)
    2351            0 :             t_control%mgrid_e_cutoff = -t_control%mgrid_e_cutoff
    2352              :          END IF
    2353              :       END IF
    2354              : 
    2355              :       ! expand XC subsection (if given explicitly)
    2356         9064 :       xc_section => section_vals_get_subs_vals(t_section, "XC")
    2357         9064 :       xc_func => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
    2358         9064 :       CALL section_vals_get(xc_func, explicit=explicit)
    2359         9064 :       IF (explicit) THEN
    2360          298 :          CALL xc_functionals_expand(xc_func, xc_section)
    2361              :       END IF
    2362              : 
    2363              :       ! sTDA subsection
    2364         9064 :       stda_section => section_vals_get_subs_vals(t_section, "STDA")
    2365         9064 :       IF (t_control%kernel == tddfpt_kernel_stda) THEN
    2366          132 :          t_control%stda_control%hfx_fraction = 0.0_dp
    2367          132 :          t_control%stda_control%do_exchange = .TRUE.
    2368          132 :          t_control%stda_control%eps_td_filter = 1.e-10_dp
    2369          132 :          t_control%stda_control%mn_alpha = -99.0_dp
    2370          132 :          t_control%stda_control%mn_beta = -99.0_dp
    2371              :          ! set default for Ewald method (on/off) dependent on periodicity
    2372          236 :          SELECT CASE (qs_control%periodicity)
    2373              :          CASE (0)
    2374          104 :             t_control%stda_control%do_ewald = .FALSE.
    2375              :          CASE (1)
    2376            0 :             t_control%stda_control%do_ewald = .TRUE.
    2377              :          CASE (2)
    2378            0 :             t_control%stda_control%do_ewald = .TRUE.
    2379              :          CASE (3)
    2380           28 :             t_control%stda_control%do_ewald = .TRUE.
    2381              :          CASE DEFAULT
    2382          132 :             CPABORT("Illegal value for periodiciy")
    2383              :          END SELECT
    2384          132 :          CALL section_vals_get(stda_section, explicit=explicit)
    2385          132 :          IF (explicit) THEN
    2386          116 :             CALL section_vals_val_get(stda_section, "HFX_FRACTION", r_val=hfx, explicit=expl)
    2387          116 :             IF (expl) t_control%stda_control%hfx_fraction = hfx
    2388          116 :             CALL section_vals_val_get(stda_section, "EPS_TD_FILTER", r_val=filter, explicit=expl)
    2389          116 :             IF (expl) t_control%stda_control%eps_td_filter = filter
    2390          116 :             CALL section_vals_val_get(stda_section, "DO_EWALD", l_val=do_ewald, explicit=expl)
    2391          116 :             IF (expl) t_control%stda_control%do_ewald = do_ewald
    2392          116 :             CALL section_vals_val_get(stda_section, "DO_EXCHANGE", l_val=do_exchange, explicit=expl)
    2393          116 :             IF (expl) t_control%stda_control%do_exchange = do_exchange
    2394          116 :             CALL section_vals_val_get(stda_section, "MATAGA_NISHIMOTO_CEXP", r_val=fval)
    2395          116 :             t_control%stda_control%mn_alpha = fval
    2396          116 :             CALL section_vals_val_get(stda_section, "MATAGA_NISHIMOTO_XEXP", r_val=fval)
    2397          116 :             t_control%stda_control%mn_beta = fval
    2398              :          END IF
    2399          132 :          CALL section_vals_val_get(stda_section, "COULOMB_SR_CUT", r_val=fval)
    2400          132 :          t_control%stda_control%coulomb_sr_cut = fval
    2401          132 :          CALL section_vals_val_get(stda_section, "COULOMB_SR_EPS", r_val=fval)
    2402          132 :          t_control%stda_control%coulomb_sr_eps = fval
    2403              :       END IF
    2404              : 
    2405         9064 :       CALL timestop(handle)
    2406         9064 :    END SUBROUTINE read_tddfpt2_control
    2407              : 
    2408              : ! **************************************************************************************************
    2409              : !> \brief Write the DFT control parameters to the output unit.
    2410              : !> \param dft_control ...
    2411              : !> \param dft_section ...
    2412              : ! **************************************************************************************************
    2413        15506 :    SUBROUTINE write_dft_control(dft_control, dft_section)
    2414              :       TYPE(dft_control_type), POINTER                    :: dft_control
    2415              :       TYPE(section_vals_type), POINTER                   :: dft_section
    2416              : 
    2417              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'write_dft_control'
    2418              : 
    2419              :       CHARACTER(LEN=20)                                  :: tmpStr
    2420              :       INTEGER                                            :: handle, i, i_rep, max_mtlr_iter, n_rep, &
    2421              :                                                             output_unit
    2422              :       REAL(kind=dp)                                      :: density_cut, density_smooth_cut_range, &
    2423              :                                                             eps_u_j_loop, gradient_cut, tau_cut
    2424              :       TYPE(cp_logger_type), POINTER                      :: logger
    2425              :       TYPE(enumeration_type), POINTER                    :: enum
    2426              :       TYPE(keyword_type), POINTER                        :: keyword
    2427              :       TYPE(section_type), POINTER                        :: section
    2428              :       TYPE(section_vals_type), POINTER                   :: xc_section
    2429              : 
    2430        10536 :       IF (dft_control%qs_control%semi_empirical) RETURN
    2431         8014 :       IF (dft_control%qs_control%dftb) RETURN
    2432         7716 :       IF (dft_control%qs_control%xtb) THEN
    2433         1224 :          CALL write_xtb_control(dft_control%qs_control%xtb_control, dft_section)
    2434         1224 :          RETURN
    2435              :       END IF
    2436         6492 :       CALL timeset(routineN, handle)
    2437              : 
    2438         6492 :       NULLIFY (logger)
    2439         6492 :       logger => cp_get_default_logger()
    2440              : 
    2441              :       output_unit = cp_print_key_unit_nr(logger, dft_section, &
    2442         6492 :                                          "PRINT%DFT_CONTROL_PARAMETERS", extension=".Log")
    2443              : 
    2444         6492 :       IF (output_unit > 0) THEN
    2445              : 
    2446         1552 :          xc_section => section_vals_get_subs_vals(dft_section, "XC")
    2447              : 
    2448         1552 :          IF (dft_control%uks) THEN
    2449              :             WRITE (UNIT=output_unit, FMT="(/,T2,A,T78,A)") &
    2450          441 :                "DFT| Spin unrestricted (spin-polarized) Kohn-Sham calculation", "UKS"
    2451         1111 :          ELSE IF (dft_control%roks) THEN
    2452              :             WRITE (UNIT=output_unit, FMT="(/,T2,A,T77,A)") &
    2453           15 :                "DFT| Spin restricted open Kohn-Sham calculation", "ROKS"
    2454              :          ELSE
    2455              :             WRITE (UNIT=output_unit, FMT="(/,T2,A,T78,A)") &
    2456         1096 :                "DFT| Spin restricted Kohn-Sham (RKS) calculation", "RKS"
    2457              :          END IF
    2458              : 
    2459              :          WRITE (UNIT=output_unit, FMT="(T2,A,T76,I5)") &
    2460         1552 :             "DFT| Multiplicity", dft_control%multiplicity
    2461              :          WRITE (UNIT=output_unit, FMT="(T2,A,T76,I5)") &
    2462         1552 :             "DFT| Number of spin states", dft_control%nspins
    2463              : 
    2464              :          WRITE (UNIT=output_unit, FMT="(T2,A,T76,I5)") &
    2465         1552 :             "DFT| Charge", dft_control%charge
    2466              : 
    2467         1552 :          IF (dft_control%sic_method_id /= sic_none) CALL cite_reference(VandeVondele2005b)
    2468         3090 :          SELECT CASE (dft_control%sic_method_id)
    2469              :          CASE (sic_none)
    2470         1538 :             tmpstr = "NO"
    2471              :          CASE (sic_mauri_spz)
    2472            6 :             tmpstr = "SPZ/MAURI SIC"
    2473              :          CASE (sic_mauri_us)
    2474            3 :             tmpstr = "US/MAURI SIC"
    2475              :          CASE (sic_ad)
    2476            3 :             tmpstr = "AD SIC"
    2477              :          CASE (sic_eo)
    2478            2 :             tmpstr = "Explicit Orbital SIC"
    2479              :          CASE DEFAULT
    2480              :             ! fix throughout the cp2k for this option
    2481         1552 :             CPABORT("SIC option unknown")
    2482              :          END SELECT
    2483              : 
    2484              :          WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") &
    2485         1552 :             "DFT| Self-interaction correction (SIC)", ADJUSTR(TRIM(tmpstr))
    2486              : 
    2487         1552 :          IF (dft_control%sic_method_id /= sic_none) THEN
    2488              :             WRITE (UNIT=output_unit, FMT="(T2,A,T66,ES15.6)") &
    2489           14 :                "DFT| SIC scaling parameter a", dft_control%sic_scaling_a, &
    2490           28 :                "DFT| SIC scaling parameter b", dft_control%sic_scaling_b
    2491              :          END IF
    2492              : 
    2493         1552 :          IF (dft_control%sic_method_id == sic_eo) THEN
    2494            2 :             IF (dft_control%sic_list_id == sic_list_all) THEN
    2495              :                WRITE (UNIT=output_unit, FMT="(T2,A,T66,A)") &
    2496            1 :                   "DFT| SIC orbitals", "ALL"
    2497              :             END IF
    2498            2 :             IF (dft_control%sic_list_id == sic_list_unpaired) THEN
    2499              :                WRITE (UNIT=output_unit, FMT="(T2,A,T66,A)") &
    2500            1 :                   "DFT| SIC orbitals", "UNPAIRED"
    2501              :             END IF
    2502              :          END IF
    2503              : 
    2504         1552 :          CALL section_vals_val_get(xc_section, "density_cutoff", r_val=density_cut)
    2505         1552 :          CALL section_vals_val_get(xc_section, "gradient_cutoff", r_val=gradient_cut)
    2506         1552 :          CALL section_vals_val_get(xc_section, "tau_cutoff", r_val=tau_cut)
    2507         1552 :          CALL section_vals_val_get(xc_section, "density_smooth_cutoff_range", r_val=density_smooth_cut_range)
    2508              : 
    2509              :          WRITE (UNIT=output_unit, FMT="(T2,A,T66,ES15.6)") &
    2510         1552 :             "DFT| Cutoffs: density ", density_cut, &
    2511         1552 :             "DFT|          gradient", gradient_cut, &
    2512         1552 :             "DFT|          tau     ", tau_cut, &
    2513         3104 :             "DFT|          cutoff_smoothing_range", density_smooth_cut_range
    2514              :          CALL section_vals_val_get(xc_section, "XC_GRID%XC_SMOOTH_RHO", &
    2515         1552 :                                    c_val=tmpStr)
    2516              :          WRITE (output_unit, '( A, T61, A )') &
    2517         1552 :             " DFT| XC density smoothing ", ADJUSTR(tmpStr)
    2518              :          CALL section_vals_val_get(xc_section, "XC_GRID%XC_DERIV", &
    2519         1552 :                                    c_val=tmpStr)
    2520              :          WRITE (output_unit, '( A, T61, A )') &
    2521         1552 :             " DFT| XC derivatives ", ADJUSTR(tmpStr)
    2522         1552 :          IF (dft_control%dft_plus_u) THEN
    2523           19 :             NULLIFY (enum, keyword, section)
    2524           19 :             CALL create_dft_section(section)
    2525           19 :             keyword => section_get_keyword(section, "PLUS_U_METHOD")
    2526           19 :             CALL keyword_get(keyword, enum=enum)
    2527              :             WRITE (UNIT=output_unit, FMT="(/,T2,A,T41,A40)") &
    2528           19 :                "DFT+U| Method", ADJUSTR(TRIM(enum_i2c(enum, dft_control%plus_u_method_id)))
    2529              :             WRITE (UNIT=output_unit, FMT="(T2,A)") &
    2530           19 :                "DFT+U| Check atomic kind information for details"
    2531           19 :             IF (dft_control%mtlr_u_j) THEN
    2532            0 :                CALL section_vals_val_get(dft_section, "EPS_U_J_LOOP", r_val=eps_u_j_loop)
    2533              :                WRITE (UNIT=output_unit, FMT="(T2,A,T67,ES14.7E3)") &
    2534            0 :                   "MTLR U J| EPS_U_J_LOOP", eps_u_j_loop
    2535            0 :                CALL section_vals_val_get(dft_section, "MAX_MTLR_LOOP", i_val=max_mtlr_iter)
    2536              :                WRITE (UNIT=output_unit, FMT="(T2,A,T67,I10)") &
    2537            0 :                   "MTLR U J| MAX_MTLR_LOOP", max_mtlr_iter
    2538              :             END IF
    2539           19 :             CALL section_release(section)
    2540              :          END IF
    2541              : 
    2542         1552 :          WRITE (UNIT=output_unit, FMT="(A)") ""
    2543         1552 :          CALL xc_write(output_unit, xc_section, dft_control%lsd)
    2544              : 
    2545         1552 :          IF (dft_control%apply_period_efield) THEN
    2546            6 :             WRITE (UNIT=output_unit, FMT="(A)") ""
    2547            6 :             IF (dft_control%period_efield%displacement_field) THEN
    2548              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    2549            0 :                   "PERIODIC_EFIELD| Use displacement field formulation"
    2550              :                WRITE (UNIT=output_unit, FMT="(T2,A,T66,1X,ES14.6)") &
    2551            0 :                   "PERIODIC_EFIELD| Displacement field filter: x", &
    2552            0 :                   dft_control%period_efield%d_filter(1), &
    2553            0 :                   "PERIODIC_EFIELD|                            y", &
    2554            0 :                   dft_control%period_efield%d_filter(2), &
    2555            0 :                   "PERIODIC_EFIELD|                            z", &
    2556            0 :                   dft_control%period_efield%d_filter(3)
    2557              :             END IF
    2558              :             WRITE (UNIT=output_unit, FMT="(T2,A,T66,1X,ES14.6)") &
    2559            6 :                "PERIODIC_EFIELD| Polarisation vector:       x", &
    2560            6 :                dft_control%period_efield%polarisation(1), &
    2561            6 :                "PERIODIC_EFIELD|                            y", &
    2562            6 :                dft_control%period_efield%polarisation(2), &
    2563            6 :                "PERIODIC_EFIELD|                            z", &
    2564           12 :                dft_control%period_efield%polarisation(3)
    2565              : 
    2566              :             WRITE (UNIT=output_unit, FMT="(T2,A,T66,1X,I14)") &
    2567            6 :                "PERIODIC_EFIELD| Start Frame:", &
    2568            6 :                dft_control%period_efield%start_frame, &
    2569            6 :                "PERIODIC_EFIELD| End Frame:", &
    2570           12 :                dft_control%period_efield%end_frame
    2571              : 
    2572            6 :             IF (ALLOCATED(dft_control%period_efield%strength_list)) THEN
    2573              :                WRITE (UNIT=output_unit, FMT="(T2,A,T66,1X,I14)") &
    2574            2 :                   "PERIODIC_EFIELD| Number of Intensities:", &
    2575            4 :                   SIZE(dft_control%period_efield%strength_list)
    2576              :                WRITE (UNIT=output_unit, FMT="(T2,A,I10,T66,1X,ES14.6)") &
    2577            2 :                   "PERIODIC_EFIELD| Intensity List [a.u.] ", &
    2578            4 :                   1, dft_control%period_efield%strength_list(1)
    2579           24 :                DO i = 2, SIZE(dft_control%period_efield%strength_list)
    2580              :                   WRITE (UNIT=output_unit, FMT="(T2,A,I10,T66,1X,ES14.6)") &
    2581           22 :                      "PERIODIC_EFIELD|                       ", &
    2582           46 :                      i, dft_control%period_efield%strength_list(i)
    2583              :                END DO
    2584              :             ELSE
    2585              :                WRITE (UNIT=output_unit, FMT="(T2,A,T66,1X,ES14.6)") &
    2586            4 :                   "PERIODIC_EFIELD| Intensity [a.u.]:", &
    2587            8 :                   dft_control%period_efield%strength
    2588              :             END IF
    2589              : 
    2590           24 :             IF (NORM2(dft_control%period_efield%polarisation) < EPSILON(0.0_dp)) THEN
    2591            0 :                CPABORT("Invalid (too small) polarisation vector specified for PERIODIC_EFIELD")
    2592              :             END IF
    2593              :          END IF
    2594              : 
    2595         1552 :          IF (dft_control%do_sccs) THEN
    2596              :             WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
    2597            5 :                "SCCS| Self-consistent continuum solvation model"
    2598              :             WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES20.6)") &
    2599            5 :                "SCCS| Relative permittivity of the solvent (medium)", &
    2600            5 :                dft_control%sccs_control%epsilon_solvent, &
    2601            5 :                "SCCS| Absolute permittivity [a.u.]", &
    2602           10 :                dft_control%sccs_control%epsilon_solvent/fourpi
    2603            9 :             SELECT CASE (dft_control%sccs_control%method_id)
    2604              :             CASE (sccs_andreussi)
    2605              :                WRITE (UNIT=output_unit, FMT="(T2,A,/,(T2,A,T61,ES20.6))") &
    2606            4 :                   "SCCS| Dielectric function proposed by Andreussi et al.", &
    2607            4 :                   "SCCS|  rho_max", dft_control%sccs_control%rho_max, &
    2608            8 :                   "SCCS|  rho_min", dft_control%sccs_control%rho_min
    2609              :             CASE (sccs_fattebert_gygi)
    2610              :                WRITE (UNIT=output_unit, FMT="(T2,A,/,(T2,A,T61,ES20.6))") &
    2611            1 :                   "SCCS| Dielectric function proposed by Fattebert and Gygi", &
    2612            1 :                   "SCCS|  beta", dft_control%sccs_control%beta, &
    2613            2 :                   "SCCS|  rho_zero", dft_control%sccs_control%rho_zero
    2614              :             CASE (sccs_saa_andreussi)
    2615              :                WRITE (UNIT=output_unit, FMT="(T2,A,/,A,/,(T2,A,T61,ES20.6))") &
    2616            0 :                   "SCCS| Dielectric function of the solvent aware algorithm", &
    2617            0 :                   "SCCS| proposed by Andreussi et al.", &
    2618            0 :                   "SCCS|  rho_max", dft_control%sccs_control%rho_max, &
    2619            0 :                   "SCCS|  rho_min", dft_control%sccs_control%rho_min, &
    2620            0 :                   "SCCS|  f0", dft_control%sccs_control%f0, &
    2621            0 :                   "SCCS|  delta_eta", dft_control%sccs_control%delta_eta, &
    2622            0 :                   "SCCS|  alpha_zeta", dft_control%sccs_control%alpha_zeta, &
    2623            0 :                   "SCCS|  delta_zeta", dft_control%sccs_control%delta_zeta, &
    2624            0 :                   "SCCS|  R_solv", dft_control%sccs_control%R_solv
    2625              :             CASE DEFAULT
    2626            5 :                CPABORT("Invalid SCCS model specified. Please, check your input!")
    2627              :             END SELECT
    2628            6 :             SELECT CASE (dft_control%sccs_control%derivative_method)
    2629              :             CASE (sccs_derivative_fft)
    2630              :                WRITE (UNIT=output_unit, FMT="(T2,A,T46,A35)") &
    2631            1 :                   "SCCS| Numerical derivative calculation", &
    2632            2 :                   ADJUSTR("FFT")
    2633              :             CASE (sccs_derivative_cd3)
    2634              :                WRITE (UNIT=output_unit, FMT="(T2,A,T46,A35)") &
    2635            0 :                   "SCCS| Numerical derivative calculation", &
    2636            0 :                   ADJUSTR("3-point stencil central differences")
    2637              :             CASE (sccs_derivative_cd5)
    2638              :                WRITE (UNIT=output_unit, FMT="(T2,A,T46,A35)") &
    2639            4 :                   "SCCS| Numerical derivative calculation", &
    2640            8 :                   ADJUSTR("5-point stencil central differences")
    2641              :             CASE (sccs_derivative_cd7)
    2642              :                WRITE (UNIT=output_unit, FMT="(T2,A,T46,A35)") &
    2643            0 :                   "SCCS| Numerical derivative calculation", &
    2644            0 :                   ADJUSTR("7-point stencil central differences")
    2645              :             CASE DEFAULT
    2646              :                CALL cp_abort(__LOCATION__, &
    2647              :                              "Invalid derivative method specified for SCCS model. "// &
    2648            5 :                              "Please, check your input!")
    2649              :             END SELECT
    2650              :             WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES20.6)") &
    2651            5 :                "SCCS| Repulsion parameter alpha [mN/m] = [dyn/cm]", &
    2652           10 :                cp_unit_from_cp2k(dft_control%sccs_control%alpha_solvent, "mN/m")
    2653              :             WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES20.6)") &
    2654            5 :                "SCCS| Dispersion parameter beta [GPa]", &
    2655           10 :                cp_unit_from_cp2k(dft_control%sccs_control%beta_solvent, "GPa")
    2656              :             WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES20.6)") &
    2657            5 :                "SCCS| Surface tension gamma [mN/m] = [dyn/cm]", &
    2658           10 :                cp_unit_from_cp2k(dft_control%sccs_control%gamma_solvent, "mN/m")
    2659              :             WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES20.6)") &
    2660            5 :                "SCCS| Mixing parameter applied during the iteration cycle", &
    2661           10 :                dft_control%sccs_control%mixing
    2662              :             WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES20.6)") &
    2663            5 :                "SCCS| Tolerance for the convergence of the SCCS iteration cycle", &
    2664           10 :                dft_control%sccs_control%eps_sccs
    2665              :             WRITE (UNIT=output_unit, FMT="(T2,A,T61,I20)") &
    2666            5 :                "SCCS| Maximum number of iteration steps", &
    2667           10 :                dft_control%sccs_control%max_iter
    2668              :             WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES20.6)") &
    2669            5 :                "SCCS| SCF convergence threshold for starting the SCCS iteration", &
    2670           10 :                dft_control%sccs_control%eps_scf
    2671              :             WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES20.6)") &
    2672            5 :                "SCCS| Numerical increment for the cavity surface calculation", &
    2673           10 :                dft_control%sccs_control%delta_rho
    2674              :          END IF
    2675              : 
    2676         1552 :          WRITE (UNIT=output_unit, FMT="(A)") ""
    2677              : 
    2678              :       END IF
    2679              : 
    2680         6492 :       IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
    2681            4 :          n_rep = SIZE(dft_control%probe)
    2682            4 :          IF (output_unit > 0) THEN
    2683            6 :             DO i_rep = 1, n_rep
    2684              :                WRITE (UNIT=output_unit, FMT="(T2,A,I5)") &
    2685            4 :                   "HP | hair probe set", i_rep
    2686              :                WRITE (UNIT=output_unit, FMT="(T2,A,T61,*(I5))") &
    2687            4 :                   "HP| atom indexes", &
    2688           12 :                   (dft_control%probe(i_rep)%atom_ids(i), i=1, dft_control%probe(i_rep)%natoms)
    2689              :                WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES20.6)") &
    2690            4 :                   "HP| potential", dft_control%probe(i_rep)%mu
    2691              :                WRITE (UNIT=output_unit, FMT="(T2,A,T61,F20.2)") &
    2692            4 :                   "HP| temperature", dft_control%probe(i_rep)%T
    2693              :                WRITE (UNIT=output_unit, FMT="(T2,A,T61,ES20.6)") &
    2694            6 :                   "HP| eps_hp", dft_control%probe(i_rep)%eps_hp
    2695              :             END DO
    2696              :          END IF
    2697              :       END IF
    2698              : 
    2699              :       CALL cp_print_key_finished_output(output_unit, logger, dft_section, &
    2700         6492 :                                         "PRINT%DFT_CONTROL_PARAMETERS")
    2701              : 
    2702         6492 :       CALL timestop(handle)
    2703              : 
    2704              :    END SUBROUTINE write_dft_control
    2705              : 
    2706              : ! **************************************************************************************************
    2707              : !> \brief Write the ADMM control parameters to the output unit.
    2708              : !> \param admm_control ...
    2709              : !> \param dft_section ...
    2710              : ! **************************************************************************************************
    2711          524 :    SUBROUTINE write_admm_control(admm_control, dft_section)
    2712              :       TYPE(admm_control_type), POINTER                   :: admm_control
    2713              :       TYPE(section_vals_type), POINTER                   :: dft_section
    2714              : 
    2715              :       INTEGER                                            :: iounit
    2716              :       TYPE(cp_logger_type), POINTER                      :: logger
    2717              : 
    2718          524 :       NULLIFY (logger)
    2719          524 :       logger => cp_get_default_logger()
    2720              : 
    2721              :       iounit = cp_print_key_unit_nr(logger, dft_section, &
    2722          524 :                                     "PRINT%DFT_CONTROL_PARAMETERS", extension=".Log")
    2723              : 
    2724          524 :       IF (iounit > 0) THEN
    2725              : 
    2726          257 :          SELECT CASE (admm_control%admm_type)
    2727              :          CASE (no_admm_type)
    2728          124 :             WRITE (UNIT=iounit, FMT="(/,T2,A,T77,A)") "ADMM| Specific ADMM type specified", "NONE"
    2729              :          CASE (admm1_type)
    2730            2 :             WRITE (UNIT=iounit, FMT="(/,T2,A,T76,A)") "ADMM| Specific ADMM type specified", "ADMM1"
    2731              :          CASE (admm2_type)
    2732            1 :             WRITE (UNIT=iounit, FMT="(/,T2,A,T76,A)") "ADMM| Specific ADMM type specified", "ADMM2"
    2733              :          CASE (admms_type)
    2734            4 :             WRITE (UNIT=iounit, FMT="(/,T2,A,T76,A)") "ADMM| Specific ADMM type specified", "ADMMS"
    2735              :          CASE (admmp_type)
    2736            1 :             WRITE (UNIT=iounit, FMT="(/,T2,A,T76,A)") "ADMM| Specific ADMM type specified", "ADMMP"
    2737              :          CASE (admmq_type)
    2738            1 :             WRITE (UNIT=iounit, FMT="(/,T2,A,T76,A)") "ADMM| Specific ADMM type specified", "ADMMQ"
    2739              :          CASE DEFAULT
    2740          133 :             CPABORT("admm_type")
    2741              :          END SELECT
    2742              : 
    2743          217 :          SELECT CASE (admm_control%purification_method)
    2744              :          CASE (do_admm_purify_none)
    2745           84 :             WRITE (UNIT=iounit, FMT="(T2,A,T77,A)") "ADMM| Density matrix purification method", "NONE"
    2746              :          CASE (do_admm_purify_cauchy)
    2747            9 :             WRITE (UNIT=iounit, FMT="(T2,A,T75,A)") "ADMM| Density matrix purification method", "Cauchy"
    2748              :          CASE (do_admm_purify_cauchy_subspace)
    2749            5 :             WRITE (UNIT=iounit, FMT="(T2,A,T66,A)") "ADMM| Density matrix purification method", "Cauchy subspace"
    2750              :          CASE (do_admm_purify_mo_diag)
    2751           25 :             WRITE (UNIT=iounit, FMT="(T2,A,T63,A)") "ADMM| Density matrix purification method", "MO diagonalization"
    2752              :          CASE (do_admm_purify_mo_no_diag)
    2753            3 :             WRITE (UNIT=iounit, FMT="(T2,A,T71,A)") "ADMM| Density matrix purification method", "MO no diag"
    2754              :          CASE (do_admm_purify_mcweeny)
    2755            1 :             WRITE (UNIT=iounit, FMT="(T2,A,T74,A)") "ADMM| Density matrix purification method", "McWeeny"
    2756              :          CASE (do_admm_purify_none_dm)
    2757            6 :             WRITE (UNIT=iounit, FMT="(T2,A,T73,A)") "ADMM| Density matrix purification method", "NONE(DM)"
    2758              :          CASE DEFAULT
    2759          133 :             CPABORT("admm_purification_method")
    2760              :          END SELECT
    2761              : 
    2762          238 :          SELECT CASE (admm_control%method)
    2763              :          CASE (do_admm_basis_projection)
    2764          105 :             WRITE (UNIT=iounit, FMT="(T2,A)") "ADMM| Orbital projection on ADMM basis"
    2765              :          CASE (do_admm_blocking_purify_full)
    2766            3 :             WRITE (UNIT=iounit, FMT="(T2,A)") "ADMM| Blocked Fock matrix projection with full purification"
    2767              :          CASE (do_admm_blocked_projection)
    2768            6 :             WRITE (UNIT=iounit, FMT="(T2,A)") "ADMM| Blocked Fock matrix projection"
    2769              :          CASE (do_admm_charge_constrained_projection)
    2770           19 :             WRITE (UNIT=iounit, FMT="(T2,A)") "ADMM| Orbital projection with charge constrain"
    2771              :          CASE DEFAULT
    2772          133 :             CPABORT("admm method")
    2773              :          END SELECT
    2774              : 
    2775          154 :          SELECT CASE (admm_control%scaling_model)
    2776              :          CASE (do_admm_exch_scaling_none)
    2777              :          CASE (do_admm_exch_scaling_merlot)
    2778           21 :             WRITE (UNIT=iounit, FMT="(T2,A)") "ADMM| Use Merlot (2014) scaling model"
    2779              :          CASE DEFAULT
    2780          133 :             CPABORT("admm scaling_model")
    2781              :          END SELECT
    2782              : 
    2783          133 :          WRITE (UNIT=iounit, FMT="(T2,A,T61,G20.10)") "ADMM| eps_filter", admm_control%eps_filter
    2784              : 
    2785          145 :          SELECT CASE (admm_control%aux_exch_func)
    2786              :          CASE (do_admm_aux_exch_func_none)
    2787           12 :             WRITE (UNIT=iounit, FMT="(T2,A)") "ADMM| No exchange functional correction term used"
    2788              :          CASE (do_admm_aux_exch_func_default, do_admm_aux_exch_func_default_libxc)
    2789           90 :             WRITE (UNIT=iounit, FMT="(T2,A,T74,A)") "ADMM| Exchange functional in correction term", "(W)PBEX"
    2790              :          CASE (do_admm_aux_exch_func_pbex, do_admm_aux_exch_func_pbex_libxc)
    2791           22 :             WRITE (UNIT=iounit, FMT="(T2,A,T77,A)") "ADMM| Exchange functional in correction term", "PBEX"
    2792              :          CASE (do_admm_aux_exch_func_opt, do_admm_aux_exch_func_opt_libxc)
    2793            8 :             WRITE (UNIT=iounit, FMT="(T2,A,T77,A)") "ADMM| Exchange functional in correction term", "OPTX"
    2794              :          CASE (do_admm_aux_exch_func_bee, do_admm_aux_exch_func_bee_libxc)
    2795            1 :             WRITE (UNIT=iounit, FMT="(T2,A,T74,A)") "ADMM| Exchange functional in correction term", "Becke88"
    2796              :          CASE (do_admm_aux_exch_func_sx_libxc)
    2797            0 :             WRITE (UNIT=iounit, FMT="(T2,A,T74,A)") "ADMM| Exchange functional in correction term", "SlaterX"
    2798              :          CASE DEFAULT
    2799          133 :             CPABORT("admm aux_exch_func")
    2800              :          END SELECT
    2801              : 
    2802          133 :          WRITE (UNIT=iounit, FMT="(A)") ""
    2803              : 
    2804              :       END IF
    2805              : 
    2806              :       CALL cp_print_key_finished_output(iounit, logger, dft_section, &
    2807          524 :                                         "PRINT%DFT_CONTROL_PARAMETERS")
    2808          524 :    END SUBROUTINE write_admm_control
    2809              : 
    2810              : ! **************************************************************************************************
    2811              : !> \brief Write the xTB control parameters to the output unit.
    2812              : !> \param xtb_control ...
    2813              : !> \param dft_section ...
    2814              : ! **************************************************************************************************
    2815         1224 :    SUBROUTINE write_xtb_control(xtb_control, dft_section)
    2816              :       TYPE(xtb_control_type), POINTER                    :: xtb_control
    2817              :       TYPE(section_vals_type), POINTER                   :: dft_section
    2818              : 
    2819              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'write_xtb_control'
    2820              : 
    2821              :       CHARACTER(LEN=16)                                  :: scc_mixer_name, solver_name
    2822              :       INTEGER                                            :: handle, output_unit
    2823              :       TYPE(cp_logger_type), POINTER                      :: logger
    2824              : 
    2825         1224 :       CALL timeset(routineN, handle)
    2826         1224 :       NULLIFY (logger)
    2827         1224 :       logger => cp_get_default_logger()
    2828              : 
    2829              :       output_unit = cp_print_key_unit_nr(logger, dft_section, &
    2830         1224 :                                          "PRINT%DFT_CONTROL_PARAMETERS", extension=".Log")
    2831              : 
    2832         1224 :       IF (output_unit > 0) THEN
    2833              : 
    2834              :          WRITE (UNIT=output_unit, FMT="(/,T2,A,T31,A50)") &
    2835          118 :             "xTB| Parameter file", ADJUSTR(TRIM(xtb_control%parameter_file_name))
    2836              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") &
    2837          118 :             "xTB| Basis expansion STO-NG", xtb_control%sto_ng
    2838              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") &
    2839          118 :             "xTB| Basis expansion STO-NG for Hydrogen", xtb_control%h_sto_ng
    2840              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,E10.4)") &
    2841          118 :             "xTB| Repulsive pair potential accuracy", xtb_control%eps_pair
    2842              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.6)") &
    2843          118 :             "xTB| Repulsive enhancement factor", xtb_control%enscale
    2844              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,L10)") &
    2845          118 :             "xTB| Halogen interaction potential", xtb_control%xb_interaction
    2846              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.3)") &
    2847          118 :             "xTB| Halogen interaction potential cutoff radius", xtb_control%xb_radius
    2848              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,L10)") &
    2849          118 :             "xTB| Nonbonded interactions", xtb_control%do_nonbonded
    2850          118 :          SELECT CASE (xtb_control%vdw_type)
    2851              :          CASE (xtb_vdw_type_none)
    2852            0 :             WRITE (UNIT=output_unit, FMT="(T2,A)") "xTB| No vdW potential selected"
    2853              :          CASE (xtb_vdw_type_d3)
    2854          117 :             WRITE (UNIT=output_unit, FMT="(T2,A,T72,A)") "xTB| vdW potential type:", "DFTD3(BJ)"
    2855              :             WRITE (UNIT=output_unit, FMT="(T2,A,T31,A50)") &
    2856          117 :                "xTB| D3 Dispersion: Parameter file", ADJUSTR(TRIM(xtb_control%dispersion_parameter_file))
    2857              :          CASE (xtb_vdw_type_d4)
    2858            1 :             WRITE (UNIT=output_unit, FMT="(T2,A,T76,A)") "xTB| vdW potential type:", "DFTD4"
    2859              :             WRITE (UNIT=output_unit, FMT="(T2,A,T31,A50)") &
    2860            1 :                "xTB| D4 Dispersion: Parameter file", ADJUSTR(TRIM(xtb_control%dispersion_parameter_file))
    2861              :          CASE DEFAULT
    2862          118 :             CPABORT("vdw type")
    2863              :          END SELECT
    2864              :          WRITE (UNIT=output_unit, FMT="(T2,A,T51,3F10.3)") &
    2865          118 :             "xTB| Huckel constants ks kp kd", xtb_control%ks, xtb_control%kp, xtb_control%kd
    2866              :          WRITE (UNIT=output_unit, FMT="(T2,A,T61,2F10.3)") &
    2867          118 :             "xTB| Huckel constants ksp k2sh", xtb_control%ksp, xtb_control%k2sh
    2868              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.3)") &
    2869          118 :             "xTB| Mataga-Nishimoto exponent", xtb_control%kg
    2870              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.3)") &
    2871          118 :             "xTB| Repulsion potential exponent", xtb_control%kf
    2872              :          WRITE (UNIT=output_unit, FMT="(T2,A,T51,3F10.3)") &
    2873          118 :             "xTB| Coordination number scaling kcn(s) kcn(p) kcn(d)", &
    2874          236 :             xtb_control%kcns, xtb_control%kcnp, xtb_control%kcnd
    2875              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.3)") &
    2876          118 :             "xTB| Electronegativity scaling", xtb_control%ken
    2877              :          WRITE (UNIT=output_unit, FMT="(T2,A,T61,2F10.3)") &
    2878          118 :             "xTB| Halogen potential scaling kxr kx2", xtb_control%kxr, xtb_control%kx2
    2879          207 :          SELECT CASE (xtb_control%tblite_scc_mixer)
    2880              :          CASE (tblite_scc_mixer_auto)
    2881           89 :             scc_mixer_name = "AUTO"
    2882              :          CASE (tblite_scc_mixer_tblite)
    2883            8 :             scc_mixer_name = "TBLITE"
    2884              :          CASE (tblite_scc_mixer_cp2k)
    2885            4 :             scc_mixer_name = "CP2K"
    2886              :          CASE (tblite_scc_mixer_none)
    2887           17 :             scc_mixer_name = "NONE"
    2888              :          CASE DEFAULT
    2889          118 :             CPABORT("Unknown tblite SCC mixer")
    2890              :          END SELECT
    2891          236 :          SELECT CASE (xtb_control%tblite_mixer_solver)
    2892              :          CASE (tblite_solver_gvd)
    2893          118 :             solver_name = "GVD"
    2894              :          CASE (tblite_solver_gvr)
    2895            0 :             solver_name = "GVR"
    2896              :          CASE DEFAULT
    2897          118 :             CPABORT("Unknown tblite SCC mixer solver")
    2898              :          END SELECT
    2899              :          WRITE (UNIT=output_unit, FMT="(T2,A,T72,A)") &
    2900          118 :             "xTB| SCC mixer:", TRIM(scc_mixer_name)
    2901              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,ES10.3)") &
    2902          118 :             "xTB| tblite accuracy:", xtb_control%tblite_accuracy
    2903          118 :          IF (LEN_TRIM(xtb_control%tblite_param_file) > 0) THEN
    2904              :             WRITE (UNIT=output_unit, FMT="(T2,A,T33,A)") &
    2905            0 :                "xTB| tblite parameter file:", TRIM(xtb_control%tblite_param_file)
    2906              :          END IF
    2907              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.3)") &
    2908          118 :             "xTB| tblite SCC mixer damping:", xtb_control%tblite_mixer_damping
    2909              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") &
    2910          118 :             "xTB| tblite SCC mixer iterations:", xtb_control%tblite_mixer_iterations
    2911              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") &
    2912          118 :             "xTB| tblite SCC mixer memory:", xtb_control%tblite_mixer_memory
    2913              :          WRITE (UNIT=output_unit, FMT="(T2,A,T72,A)") &
    2914          118 :             "xTB| tblite SCC mixer solver:", TRIM(solver_name)
    2915              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,ES10.3)") &
    2916          118 :             "xTB| tblite SCC mixer omega0:", xtb_control%tblite_mixer_omega0
    2917              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,ES10.3)") &
    2918          118 :             "xTB| tblite SCC mixer min weight:", xtb_control%tblite_mixer_min_weight
    2919              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,ES10.3)") &
    2920          118 :             "xTB| tblite SCC mixer max weight:", xtb_control%tblite_mixer_max_weight
    2921              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,ES10.3)") &
    2922          118 :             "xTB| tblite SCC mixer weight factor:", xtb_control%tblite_mixer_weight_factor
    2923          118 :          WRITE (UNIT=output_unit, FMT="(/)")
    2924              : 
    2925              :       END IF
    2926              : 
    2927              :       CALL cp_print_key_finished_output(output_unit, logger, dft_section, &
    2928         1224 :                                         "PRINT%DFT_CONTROL_PARAMETERS")
    2929              : 
    2930         1224 :       CALL timestop(handle)
    2931              : 
    2932         1224 :    END SUBROUTINE write_xtb_control
    2933              : 
    2934              : ! **************************************************************************************************
    2935              : !> \brief Purpose: Write the QS control parameters to the output unit.
    2936              : !> \param qs_control ...
    2937              : !> \param dft_section ...
    2938              : ! **************************************************************************************************
    2939        15506 :    SUBROUTINE write_qs_control(qs_control, dft_section)
    2940              :       TYPE(qs_control_type), INTENT(IN)                  :: qs_control
    2941              :       TYPE(section_vals_type), POINTER                   :: dft_section
    2942              : 
    2943              :       CHARACTER(len=*), PARAMETER                        :: routineN = 'write_qs_control'
    2944              : 
    2945              :       CHARACTER(len=20)                                  :: method, quadrature
    2946              :       INTEGER                                            :: handle, i, igrid_level, ngrid_level, &
    2947              :                                                             output_unit
    2948              :       TYPE(cp_logger_type), POINTER                      :: logger
    2949              :       TYPE(ddapc_restraint_type), POINTER                :: ddapc_restraint_control
    2950              :       TYPE(enumeration_type), POINTER                    :: enum
    2951              :       TYPE(keyword_type), POINTER                        :: keyword
    2952              :       TYPE(section_type), POINTER                        :: qs_section
    2953              :       TYPE(section_vals_type), POINTER                   :: print_section_vals, qs_section_vals
    2954              : 
    2955        10536 :       IF (qs_control%semi_empirical) RETURN
    2956         8014 :       IF (qs_control%dftb) RETURN
    2957         7716 :       IF (qs_control%xtb) RETURN
    2958         6492 :       CALL timeset(routineN, handle)
    2959         6492 :       NULLIFY (logger, print_section_vals, qs_section, qs_section_vals)
    2960         6492 :       logger => cp_get_default_logger()
    2961         6492 :       print_section_vals => section_vals_get_subs_vals(dft_section, "PRINT")
    2962         6492 :       qs_section_vals => section_vals_get_subs_vals(dft_section, "QS")
    2963         6492 :       CALL section_vals_get(qs_section_vals, section=qs_section)
    2964              : 
    2965         6492 :       NULLIFY (enum, keyword)
    2966         6492 :       keyword => section_get_keyword(qs_section, "METHOD")
    2967         6492 :       CALL keyword_get(keyword, enum=enum)
    2968         6492 :       method = TRIM(enum_i2c(enum, qs_control%method_id))
    2969              : 
    2970         6492 :       NULLIFY (enum, keyword)
    2971         6492 :       keyword => section_get_keyword(qs_section, "QUADRATURE")
    2972         6492 :       CALL keyword_get(keyword, enum=enum)
    2973         6492 :       quadrature = TRIM(enum_i2c(enum, qs_control%gapw_control%quadrature))
    2974              : 
    2975              :       output_unit = cp_print_key_unit_nr(logger, print_section_vals, &
    2976         6492 :                                          "DFT_CONTROL_PARAMETERS", extension=".Log")
    2977         6492 :       IF (output_unit > 0) THEN
    2978         1552 :          ngrid_level = SIZE(qs_control%e_cutoff)
    2979              :          WRITE (UNIT=output_unit, FMT="(/,T2,A,T61,A20)") &
    2980         1552 :             "QS| Method:", ADJUSTR(method)
    2981         1552 :          IF (qs_control%pw_grid_opt%spherical) THEN
    2982              :             WRITE (UNIT=output_unit, FMT="(T2,A,T61,A)") &
    2983            0 :                "QS| Density plane wave grid type", " SPHERICAL HALFSPACE"
    2984         1552 :          ELSE IF (qs_control%pw_grid_opt%fullspace) THEN
    2985              :             WRITE (UNIT=output_unit, FMT="(T2,A,T57,A)") &
    2986         1552 :                "QS| Density plane wave grid type", " NON-SPHERICAL FULLSPACE"
    2987              :          ELSE
    2988              :             WRITE (UNIT=output_unit, FMT="(T2,A,T57,A)") &
    2989            0 :                "QS| Density plane wave grid type", " NON-SPHERICAL HALFSPACE"
    2990              :          END IF
    2991              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") &
    2992         1552 :             "QS| Number of grid levels:", SIZE(qs_control%e_cutoff)
    2993         1552 :          IF (ngrid_level == 1) THEN
    2994              :             WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.1)") &
    2995           80 :                "QS| Density cutoff [a.u.]:", qs_control%e_cutoff(1)
    2996              :          ELSE
    2997              :             WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.1)") &
    2998         1472 :                "QS| Density cutoff [a.u.]:", qs_control%cutoff
    2999         1472 :             IF (qs_control%commensurate_mgrids) THEN
    3000          132 :                WRITE (UNIT=output_unit, FMT="(T2,A)") "QS| Using commensurate multigrids"
    3001              :             END IF
    3002              :             WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.1)") &
    3003         1472 :                "QS| Multi grid cutoff [a.u.]: 1) grid level", qs_control%e_cutoff(1)
    3004              :             WRITE (UNIT=output_unit, FMT="(T2,A,I3,A,T71,F10.1)") &
    3005         4584 :                ("QS|                         ", igrid_level, ") grid level", &
    3006         6056 :                 qs_control%e_cutoff(igrid_level), &
    3007         7528 :                 igrid_level=2, SIZE(qs_control%e_cutoff))
    3008              :          END IF
    3009         1552 :          IF (qs_control%pao) THEN
    3010            0 :             WRITE (UNIT=output_unit, FMT="(T2,A)") "QS| PAO active"
    3011              :          END IF
    3012              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.1)") &
    3013         1552 :             "QS| Grid level progression factor:", qs_control%progression_factor
    3014              :          WRITE (UNIT=output_unit, FMT="(T2,A,T71,F10.1)") &
    3015         1552 :             "QS| Relative density cutoff [a.u.]:", qs_control%relative_cutoff
    3016              :          WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") &
    3017         1552 :             "QS| Interaction thresholds: eps_pgf_orb:", &
    3018         1552 :             qs_control%eps_pgf_orb, &
    3019         1552 :             "QS|                         eps_filter_matrix:", &
    3020         1552 :             qs_control%eps_filter_matrix, &
    3021         1552 :             "QS|                         eps_core_charge:", &
    3022         1552 :             qs_control%eps_core_charge, &
    3023         1552 :             "QS|                         eps_rho_gspace:", &
    3024         1552 :             qs_control%eps_rho_gspace, &
    3025         1552 :             "QS|                         eps_rho_rspace:", &
    3026         1552 :             qs_control%eps_rho_rspace, &
    3027         1552 :             "QS|                         eps_gvg_rspace:", &
    3028         1552 :             qs_control%eps_gvg_rspace, &
    3029         1552 :             "QS|                         eps_ppl:", &
    3030         1552 :             qs_control%eps_ppl, &
    3031         1552 :             "QS|                         eps_ppnl:", &
    3032         3104 :             qs_control%eps_ppnl
    3033         1552 :          IF (qs_control%gapw) THEN
    3034          293 :             IF (qs_control%gapw_control%accurate_xcint) THEN
    3035              :                WRITE (UNIT=output_unit, FMT="(T2,A,T69,F12.6)") &
    3036           47 :                   "QS| GAPW|      XC integration using accurate scheme: Ref. exponent =", &
    3037           94 :                   qs_control%gapw_control%aweights
    3038              :             END IF
    3039              :             !
    3040          561 :             SELECT CASE (qs_control%gapw_control%basis_1c)
    3041              :             CASE (gapw_1c_orb)
    3042              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    3043          268 :                   "QS| GAPW|      One center basis from orbital basis primitives"
    3044              :             CASE (gapw_1c_small)
    3045              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    3046           21 :                   "QS| GAPW|      One center basis extended with primitives (small:s)"
    3047              :             CASE (gapw_1c_medium)
    3048              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    3049            1 :                   "QS| GAPW|      One center basis extended with primitives (medium:sp)"
    3050              :             CASE (gapw_1c_large)
    3051              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    3052            2 :                   "QS| GAPW|      One center basis extended with primitives (large:spd)"
    3053              :             CASE (gapw_1c_very_large)
    3054              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    3055            1 :                   "QS| GAPW|      One center basis extended with primitives (very large:spdf)"
    3056              :             CASE DEFAULT
    3057          293 :                CPABORT("basis_1c incorrect")
    3058              :             END SELECT
    3059              :             WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") &
    3060          293 :                "QS| GAPW|                   eps_fit:", &
    3061          293 :                qs_control%gapw_control%eps_fit, &
    3062          293 :                "QS| GAPW|                   eps_iso:", &
    3063          293 :                qs_control%gapw_control%eps_iso, &
    3064          293 :                "QS| GAPW|                   eps_svd:", &
    3065          293 :                qs_control%gapw_control%eps_svd, &
    3066          293 :                "QS| GAPW|                   eps_cpc:", &
    3067          586 :                qs_control%gapw_control%eps_cpc
    3068              :             WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") &
    3069          293 :                "QS| GAPW|   atom-r-grid: quadrature:", &
    3070          586 :                ADJUSTR(quadrature)
    3071              :             WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") &
    3072          293 :                "QS| GAPW|      atom-s-grid:  max l :", &
    3073          293 :                qs_control%gapw_control%lmax_sphere, &
    3074          293 :                "QS| GAPW|      max_l_rho0 :", &
    3075          586 :                qs_control%gapw_control%lmax_rho0
    3076          293 :             IF (qs_control%gapw_control%non_paw_atoms) THEN
    3077              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    3078           62 :                   "QS| GAPW|      At least one kind is NOT PAW, i.e. it has only soft AO "
    3079              :             END IF
    3080          293 :             IF (qs_control%gapw_control%nopaw_as_gpw) THEN
    3081              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    3082           62 :                   "QS| GAPW|      The NOT PAW atoms are treated fully GPW"
    3083              :             END IF
    3084              :          END IF
    3085         1552 :          IF (qs_control%gapw_xc) THEN
    3086           71 :             SELECT CASE (qs_control%gapw_control%basis_1c)
    3087              :             CASE (gapw_1c_orb)
    3088              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    3089           32 :                   "QS| GAPW_XC|      One center basis from orbital basis primitives"
    3090              :             CASE (gapw_1c_small)
    3091              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    3092            7 :                   "QS| GAPW_XC|      One center basis extended with primitives (small:s)"
    3093              :             CASE (gapw_1c_medium)
    3094              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    3095            0 :                   "QS| GAPW_XC|      One center basis extended with primitives (medium:sp)"
    3096              :             CASE (gapw_1c_large)
    3097              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    3098            0 :                   "QS| GAPW_XC|      One center basis extended with primitives (large:spd)"
    3099              :             CASE (gapw_1c_very_large)
    3100              :                WRITE (UNIT=output_unit, FMT="(T2,A)") &
    3101            0 :                   "QS| GAPW_XC|      One center basis extended with primitives (very large:spdf)"
    3102              :             CASE DEFAULT
    3103           39 :                CPABORT("basis_1c incorrect")
    3104              :             END SELECT
    3105              :             WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") &
    3106           39 :                "QS| GAPW_XC|                eps_fit:", &
    3107           39 :                qs_control%gapw_control%eps_fit, &
    3108           39 :                "QS| GAPW_XC|                eps_iso:", &
    3109           39 :                qs_control%gapw_control%eps_iso, &
    3110           39 :                "QS| GAPW_XC|                eps_svd:", &
    3111           78 :                qs_control%gapw_control%eps_svd
    3112              :             WRITE (UNIT=output_unit, FMT="(T2,A,T55,A30)") &
    3113           39 :                "QS| GAPW_XC|atom-r-grid: quadrature:", &
    3114           78 :                enum_i2c(enum, qs_control%gapw_control%quadrature)
    3115              :             WRITE (UNIT=output_unit, FMT="(T2,A,T71,I10)") &
    3116           39 :                "QS| GAPW_XC|   atom-s-grid:  max l :", &
    3117           78 :                qs_control%gapw_control%lmax_sphere
    3118              :          END IF
    3119         1552 :          IF (qs_control%mulliken_restraint) THEN
    3120              :             WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") &
    3121            1 :                "QS| Mulliken restraint target", qs_control%mulliken_restraint_control%target
    3122              :             WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") &
    3123            1 :                "QS| Mulliken restraint strength", qs_control%mulliken_restraint_control%strength
    3124              :             WRITE (UNIT=output_unit, FMT="(T2,A,T73,I8)") &
    3125            1 :                "QS| Mulliken restraint atoms: ", qs_control%mulliken_restraint_control%natoms
    3126            2 :             WRITE (UNIT=output_unit, FMT="(5I8)") qs_control%mulliken_restraint_control%atoms
    3127              :          END IF
    3128         1552 :          IF (qs_control%ddapc_restraint) THEN
    3129           14 :             DO i = 1, SIZE(qs_control%ddapc_restraint_control)
    3130            8 :                ddapc_restraint_control => qs_control%ddapc_restraint_control(i)
    3131            8 :                IF (SIZE(qs_control%ddapc_restraint_control) > 1) THEN
    3132              :                   WRITE (UNIT=output_unit, FMT="(T2,A,T3,I8)") &
    3133            3 :                      "QS| parameters for DDAPC restraint number", i
    3134              :                END IF
    3135              :                WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") &
    3136            8 :                   "QS| ddapc restraint target", ddapc_restraint_control%target
    3137              :                WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") &
    3138            8 :                   "QS| ddapc restraint strength", ddapc_restraint_control%strength
    3139              :                WRITE (UNIT=output_unit, FMT="(T2,A,T73,I8)") &
    3140            8 :                   "QS| ddapc restraint atoms: ", ddapc_restraint_control%natoms
    3141           17 :                WRITE (UNIT=output_unit, FMT="(5I8)") ddapc_restraint_control%atoms
    3142            8 :                WRITE (UNIT=output_unit, FMT="(T2,A)") "Coefficients:"
    3143           17 :                WRITE (UNIT=output_unit, FMT="(5F6.2)") ddapc_restraint_control%coeff
    3144            6 :                SELECT CASE (ddapc_restraint_control%functional_form)
    3145              :                CASE (do_ddapc_restraint)
    3146              :                   WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") &
    3147            3 :                      "QS| ddapc restraint functional form :", "RESTRAINT"
    3148              :                CASE (do_ddapc_constraint)
    3149              :                   WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") &
    3150            5 :                      "QS| ddapc restraint functional form :", "CONSTRAINT"
    3151              :                CASE DEFAULT
    3152            8 :                   CPABORT("Unknown ddapc restraint")
    3153              :                END SELECT
    3154              :             END DO
    3155              :          END IF
    3156         1552 :          IF (qs_control%s2_restraint) THEN
    3157              :             WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") &
    3158            0 :                "QS| s2 restraint target", qs_control%s2_restraint_control%target
    3159              :             WRITE (UNIT=output_unit, FMT="(T2,A,T73,ES8.1)") &
    3160            0 :                "QS| s2 restraint strength", qs_control%s2_restraint_control%strength
    3161            0 :             SELECT CASE (qs_control%s2_restraint_control%functional_form)
    3162              :             CASE (do_s2_restraint)
    3163              :                WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") &
    3164            0 :                   "QS| s2 restraint functional form :", "RESTRAINT"
    3165            0 :                CPABORT("Not yet implemented")
    3166              :             CASE (do_s2_constraint)
    3167              :                WRITE (UNIT=output_unit, FMT="(T2,A,T61,A20)") &
    3168            0 :                   "QS| s2 restraint functional form :", "CONSTRAINT"
    3169              :             CASE DEFAULT
    3170            0 :                CPABORT("Unknown ddapc restraint")
    3171              :             END SELECT
    3172              :          END IF
    3173              :       END IF
    3174              :       CALL cp_print_key_finished_output(output_unit, logger, print_section_vals, &
    3175         6492 :                                         "DFT_CONTROL_PARAMETERS")
    3176              : 
    3177         6492 :       CALL timestop(handle)
    3178              : 
    3179              :    END SUBROUTINE write_qs_control
    3180              : 
    3181              : ! **************************************************************************************************
    3182              : !> \brief reads the input parameters needed for ddapc.
    3183              : !> \param qs_control ...
    3184              : !> \param qs_section ...
    3185              : !> \param ddapc_restraint_section ...
    3186              : !> \author fschiff
    3187              : !> \note
    3188              : !>      either reads DFT%QS%DDAPC_RESTRAINT or PROPERTIES%ET_coupling
    3189              : !>      if(qs_section is present the DFT part is read, if ddapc_restraint_section
    3190              : !>      is present ET_COUPLING is read. Avoid having both!!!
    3191              : ! **************************************************************************************************
    3192           14 :    SUBROUTINE read_ddapc_section(qs_control, qs_section, ddapc_restraint_section)
    3193              : 
    3194              :       TYPE(qs_control_type), INTENT(INOUT)               :: qs_control
    3195              :       TYPE(section_vals_type), OPTIONAL, POINTER         :: qs_section, ddapc_restraint_section
    3196              : 
    3197              :       INTEGER                                            :: i, j, jj, k, n_rep
    3198           14 :       INTEGER, DIMENSION(:), POINTER                     :: tmplist
    3199           14 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: rtmplist
    3200              :       TYPE(ddapc_restraint_type), POINTER                :: ddapc_restraint_control
    3201              :       TYPE(section_vals_type), POINTER                   :: ddapc_section
    3202              : 
    3203           14 :       IF (PRESENT(ddapc_restraint_section)) THEN
    3204            0 :          IF (ASSOCIATED(qs_control%ddapc_restraint_control)) THEN
    3205            0 :             IF (SIZE(qs_control%ddapc_restraint_control) >= 2) THEN
    3206            0 :                CPABORT("ET_COUPLING cannot be used in combination with a normal restraint")
    3207              :             END IF
    3208              :          ELSE
    3209            0 :             ddapc_section => ddapc_restraint_section
    3210            0 :             ALLOCATE (qs_control%ddapc_restraint_control(1))
    3211              :          END IF
    3212              :       END IF
    3213              : 
    3214           14 :       IF (PRESENT(qs_section)) THEN
    3215           14 :          NULLIFY (ddapc_section)
    3216              :          ddapc_section => section_vals_get_subs_vals(qs_section, &
    3217           14 :                                                      "DDAPC_RESTRAINT")
    3218              :       END IF
    3219              : 
    3220           32 :       DO i = 1, SIZE(qs_control%ddapc_restraint_control)
    3221              : 
    3222           18 :          CALL ddapc_control_create(qs_control%ddapc_restraint_control(i))
    3223           18 :          ddapc_restraint_control => qs_control%ddapc_restraint_control(i)
    3224              : 
    3225              :          CALL section_vals_val_get(ddapc_section, "STRENGTH", i_rep_section=i, &
    3226           18 :                                    r_val=ddapc_restraint_control%strength)
    3227              :          CALL section_vals_val_get(ddapc_section, "TARGET", i_rep_section=i, &
    3228           18 :                                    r_val=ddapc_restraint_control%target)
    3229              :          CALL section_vals_val_get(ddapc_section, "FUNCTIONAL_FORM", i_rep_section=i, &
    3230           18 :                                    i_val=ddapc_restraint_control%functional_form)
    3231              :          CALL section_vals_val_get(ddapc_section, "ATOMS", i_rep_section=i, &
    3232           18 :                                    n_rep_val=n_rep)
    3233              :          CALL section_vals_val_get(ddapc_section, "TYPE_OF_DENSITY", i_rep_section=i, &
    3234           18 :                                    i_val=ddapc_restraint_control%density_type)
    3235              : 
    3236           18 :          jj = 0
    3237           36 :          DO k = 1, n_rep
    3238              :             CALL section_vals_val_get(ddapc_section, "ATOMS", i_rep_section=i, &
    3239           18 :                                       i_rep_val=k, i_vals=tmplist)
    3240           56 :             DO j = 1, SIZE(tmplist)
    3241           38 :                jj = jj + 1
    3242              :             END DO
    3243              :          END DO
    3244           18 :          IF (jj < 1) CPABORT("Need at least 1 atom to use ddapc constraints")
    3245           18 :          ddapc_restraint_control%natoms = jj
    3246           18 :          IF (ASSOCIATED(ddapc_restraint_control%atoms)) THEN
    3247            0 :             DEALLOCATE (ddapc_restraint_control%atoms)
    3248              :          END IF
    3249           54 :          ALLOCATE (ddapc_restraint_control%atoms(ddapc_restraint_control%natoms))
    3250           18 :          jj = 0
    3251           36 :          DO k = 1, n_rep
    3252              :             CALL section_vals_val_get(ddapc_section, "ATOMS", i_rep_section=i, &
    3253           18 :                                       i_rep_val=k, i_vals=tmplist)
    3254           56 :             DO j = 1, SIZE(tmplist)
    3255           20 :                jj = jj + 1
    3256           38 :                ddapc_restraint_control%atoms(jj) = tmplist(j)
    3257              :             END DO
    3258              :          END DO
    3259              : 
    3260           18 :          IF (ASSOCIATED(ddapc_restraint_control%coeff)) THEN
    3261            0 :             DEALLOCATE (ddapc_restraint_control%coeff)
    3262              :          END IF
    3263           54 :          ALLOCATE (ddapc_restraint_control%coeff(ddapc_restraint_control%natoms))
    3264           38 :          ddapc_restraint_control%coeff = 1.0_dp
    3265              : 
    3266              :          CALL section_vals_val_get(ddapc_section, "COEFF", i_rep_section=i, &
    3267           18 :                                    n_rep_val=n_rep)
    3268           18 :          jj = 0
    3269           20 :          DO k = 1, n_rep
    3270              :             CALL section_vals_val_get(ddapc_section, "COEFF", i_rep_section=i, &
    3271            2 :                                       i_rep_val=k, r_vals=rtmplist)
    3272           22 :             DO j = 1, SIZE(rtmplist)
    3273            2 :                jj = jj + 1
    3274            2 :                IF (jj > ddapc_restraint_control%natoms) THEN
    3275            0 :                   CPABORT("Need the same number of coeff as there are atoms ")
    3276              :                END IF
    3277            4 :                ddapc_restraint_control%coeff(jj) = rtmplist(j)
    3278              :             END DO
    3279              :          END DO
    3280           68 :          IF (jj < ddapc_restraint_control%natoms .AND. jj /= 0) THEN
    3281            0 :             CPABORT("Need no or the same number of coeff as there are atoms.")
    3282              :          END IF
    3283              :       END DO
    3284           14 :       k = 0
    3285           32 :       DO i = 1, SIZE(qs_control%ddapc_restraint_control)
    3286           18 :          IF (qs_control%ddapc_restraint_control(i)%functional_form == &
    3287           24 :              do_ddapc_constraint) k = k + 1
    3288              :       END DO
    3289           14 :       IF (k == 2) CALL cp_abort(__LOCATION__, &
    3290            0 :                                 "Only a single constraint possible yet, try to use restraints instead ")
    3291              : 
    3292           14 :    END SUBROUTINE read_ddapc_section
    3293              : 
    3294              : ! **************************************************************************************************
    3295              : !> \brief ...
    3296              : !> \param dft_control ...
    3297              : !> \param efield_section ...
    3298              : !> \param cell ...
    3299              : ! **************************************************************************************************
    3300          342 :    SUBROUTINE read_efield_sections(dft_control, efield_section, cell)
    3301              :       TYPE(dft_control_type), POINTER                    :: dft_control
    3302              :       TYPE(section_vals_type), POINTER                   :: efield_section
    3303              :       TYPE(cell_type), OPTIONAL, POINTER                 :: cell
    3304              : 
    3305              :       CHARACTER(len=default_path_length)                 :: file_name
    3306              :       INTEGER                                            :: i, io, j, n, unit_nr
    3307              :       LOGICAL                                            :: amplitude_explicit, intensity_explicit
    3308          342 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tmp_vals
    3309              :       TYPE(efield_type), POINTER                         :: efield
    3310              :       TYPE(section_vals_type), POINTER                   :: tmp_section
    3311              : 
    3312          684 :       DO i = 1, SIZE(dft_control%efield_fields)
    3313          342 :          NULLIFY (dft_control%efield_fields(i)%efield)
    3314         1368 :          ALLOCATE (dft_control%efield_fields(i)%efield)
    3315          342 :          efield => dft_control%efield_fields(i)%efield
    3316          342 :          NULLIFY (efield%envelop_i_vars, efield%envelop_r_vars)
    3317              :          CALL section_vals_val_get(efield_section, "INTENSITY", i_rep_section=i, &
    3318          342 :                                    r_val=efield%strength, explicit=intensity_explicit)
    3319              :          CALL section_vals_val_get(efield_section, "AMPLITUDE", i_rep_section=i, &
    3320          342 :                                    r_val=efield%amplitude, explicit=amplitude_explicit)
    3321              : 
    3322          342 :          IF (intensity_explicit .AND. amplitude_explicit) THEN
    3323            0 :             CPABORT("Both INTENSITY and AMPLITUDE provided in EFIELD section.")
    3324              :          END IF
    3325              : 
    3326          342 :          IF (intensity_explicit) THEN
    3327           48 :             efield%amplitude = SQRT(efield%strength/(3.50944_dp*10.0_dp**16))
    3328              :          END IF
    3329              : 
    3330          342 :          IF (amplitude_explicit) THEN
    3331            0 :             efield%strength = (3.50944_dp*10.0_dp**16)*(efield%amplitude)**2
    3332              :          END IF
    3333              : 
    3334              :          CALL section_vals_val_get(efield_section, "POLARISATION", i_rep_section=i, &
    3335          342 :                                    r_vals=tmp_vals)
    3336         1026 :          ALLOCATE (efield%polarisation(SIZE(tmp_vals)))
    3337         2394 :          efield%polarisation = tmp_vals
    3338          342 :          IF (PRESENT(cell)) THEN
    3339          342 :             IF (ASSOCIATED(cell)) CALL cell_transform_input_cartesian(cell, efield%polarisation(1:3))
    3340              :          END IF
    3341              :          CALL section_vals_val_get(efield_section, "PHASE", i_rep_section=i, &
    3342          342 :                                    r_val=efield%phase_offset)
    3343              :          CALL section_vals_val_get(efield_section, "ENVELOP", i_rep_section=i, &
    3344          342 :                                    i_val=efield%envelop_id)
    3345              :          CALL section_vals_val_get(efield_section, "WAVELENGTH", i_rep_section=i, &
    3346          342 :                                    r_val=efield%wavelength)
    3347              :          CALL section_vals_val_get(efield_section, "VEC_POT_INITIAL", i_rep_section=i, &
    3348          342 :                                    r_vals=tmp_vals)
    3349         2394 :          efield%vec_pot_initial = tmp_vals
    3350          342 :          IF (PRESENT(cell)) THEN
    3351          342 :             IF (ASSOCIATED(cell)) CALL cell_transform_input_cartesian(cell, efield%vec_pot_initial(1:3))
    3352              :          END IF
    3353              : 
    3354         1026 :          IF (efield%envelop_id == constant_env) THEN
    3355          326 :             ALLOCATE (efield%envelop_i_vars(2))
    3356          326 :             tmp_section => section_vals_get_subs_vals(efield_section, "CONSTANT_ENV", i_rep_section=i)
    3357              :             CALL section_vals_val_get(tmp_section, "START_STEP", &
    3358          326 :                                       i_val=efield%envelop_i_vars(1))
    3359              :             CALL section_vals_val_get(tmp_section, "END_STEP", &
    3360          326 :                                       i_val=efield%envelop_i_vars(2))
    3361           16 :          ELSE IF (efield%envelop_id == gaussian_env) THEN
    3362           12 :             ALLOCATE (efield%envelop_r_vars(2))
    3363           12 :             tmp_section => section_vals_get_subs_vals(efield_section, "GAUSSIAN_ENV", i_rep_section=i)
    3364              :             CALL section_vals_val_get(tmp_section, "T0", &
    3365           12 :                                       r_val=efield%envelop_r_vars(1))
    3366              :             CALL section_vals_val_get(tmp_section, "SIGMA", &
    3367           12 :                                       r_val=efield%envelop_r_vars(2))
    3368            4 :          ELSE IF (efield%envelop_id == ramp_env) THEN
    3369            2 :             ALLOCATE (efield%envelop_i_vars(4))
    3370            2 :             tmp_section => section_vals_get_subs_vals(efield_section, "RAMP_ENV", i_rep_section=i)
    3371              :             CALL section_vals_val_get(tmp_section, "START_STEP_IN", &
    3372            2 :                                       i_val=efield%envelop_i_vars(1))
    3373              :             CALL section_vals_val_get(tmp_section, "END_STEP_IN", &
    3374            2 :                                       i_val=efield%envelop_i_vars(2))
    3375              :             CALL section_vals_val_get(tmp_section, "START_STEP_OUT", &
    3376            2 :                                       i_val=efield%envelop_i_vars(3))
    3377              :             CALL section_vals_val_get(tmp_section, "END_STEP_OUT", &
    3378            2 :                                       i_val=efield%envelop_i_vars(4))
    3379            2 :          ELSE IF (efield%envelop_id == custom_env) THEN
    3380            2 :             tmp_section => section_vals_get_subs_vals(efield_section, "CUSTOM_ENV", i_rep_section=i)
    3381            2 :             CALL section_vals_val_get(tmp_section, "EFIELD_FILE_NAME", c_val=file_name)
    3382            2 :             CALL open_file(file_name=TRIM(file_name), file_action="READ", file_status="OLD", unit_number=unit_nr)
    3383              :             !Determine the number of lines in file
    3384            2 :             n = 0
    3385           10 :             DO WHILE (.TRUE.)
    3386           12 :                READ (unit_nr, *, iostat=io)
    3387           12 :                IF (io /= 0) EXIT
    3388           10 :                n = n + 1
    3389              :             END DO
    3390            2 :             REWIND (unit_nr)
    3391            6 :             ALLOCATE (efield%envelop_r_vars(n + 1))
    3392              :             !Store the timestep of the list in the first entry of the r_vars
    3393            2 :             CALL section_vals_val_get(tmp_section, "TIMESTEP", r_val=efield%envelop_r_vars(1))
    3394              :             !Read the file
    3395           12 :             DO j = 2, n + 1
    3396           10 :                READ (unit_nr, *) efield%envelop_r_vars(j)
    3397           12 :                efield%envelop_r_vars(j) = cp_unit_to_cp2k(efield%envelop_r_vars(j), "volt/m")
    3398              :             END DO
    3399            2 :             CALL close_file(unit_nr)
    3400              :          END IF
    3401              :       END DO
    3402          342 :    END SUBROUTINE read_efield_sections
    3403              : 
    3404              : ! **************************************************************************************************
    3405              : !> \brief reads the input parameters needed real time propagation
    3406              : !> \param dft_control ...
    3407              : !> \param rtp_section ...
    3408              : !> \author fschiff
    3409              : ! **************************************************************************************************
    3410         2268 :    SUBROUTINE read_rtp_section(dft_control, rtp_section)
    3411              : 
    3412              :       TYPE(dft_control_type), INTENT(INOUT)              :: dft_control
    3413              :       TYPE(section_vals_type), POINTER                   :: rtp_section
    3414              : 
    3415              :       INTEGER                                            :: i, j, n_elems
    3416          324 :       INTEGER, DIMENSION(:), POINTER                     :: tmp
    3417              :       LOGICAL                                            :: is_present, linearize_bse_propagation, &
    3418              :                                                             local_moment_possible
    3419              :       TYPE(section_vals_type), POINTER                   :: proj_mo_section, subsection
    3420              : 
    3421         3888 :       ALLOCATE (dft_control%rtp_control)
    3422              :       CALL section_vals_val_get(rtp_section, "MAX_ITER", &
    3423          324 :                                 i_val=dft_control%rtp_control%max_iter)
    3424              :       CALL section_vals_val_get(rtp_section, "MAT_EXP", &
    3425          324 :                                 i_val=dft_control%rtp_control%mat_exp)
    3426              :       CALL section_vals_val_get(rtp_section, "ASPC_ORDER", &
    3427          324 :                                 i_val=dft_control%rtp_control%aspc_order)
    3428              :       CALL section_vals_val_get(rtp_section, "EXP_ACCURACY", &
    3429          324 :                                 r_val=dft_control%rtp_control%eps_exp)
    3430              :       CALL section_vals_val_get(rtp_section, "RTBSE%_SECTION_PARAMETERS_", &
    3431          324 :                                 i_val=dft_control%rtp_control%rtp_method)
    3432              :       CALL section_vals_val_get(rtp_section, "RTBSE%RTBSE_HAMILTONIAN", &
    3433          324 :                                 i_val=dft_control%rtp_control%rtbse_ham)
    3434              :       CALL section_vals_val_get(rtp_section, "RTBSE%LINEARIZED_BSE_PROPAGATION", &
    3435          324 :                                 l_val=linearize_bse_propagation)
    3436              :       ! Change rtp_method to linearized bse. The section parameter also feeds bs_env%rtp_method,
    3437              :       ! which gates the W(w=0) build in the GW step - TDDFT there would dispatch the linearized
    3438              :       ! propagator with no screened interaction to propagate with, so reject the combination.
    3439          324 :       IF (linearize_bse_propagation) THEN
    3440           58 :          IF (dft_control%rtp_control%rtp_method /= rtp_method_bse) THEN
    3441              :             CALL cp_abort(__LOCATION__, &
    3442              :                           "LINEARIZED_BSE_PROPAGATION requires the RTBSE section opened as "// &
    3443            0 :                           "'&RTBSE' or '&RTBSE RTBSE', not '&RTBSE TDDFT'.")
    3444              :          END IF
    3445           58 :          dft_control%rtp_control%rtp_method = rtp_method_bse_linearized
    3446              :       END IF
    3447              : 
    3448              :       CALL section_vals_val_get(rtp_section, "PROPAGATOR", &
    3449          324 :                                 i_val=dft_control%rtp_control%propagator)
    3450              :       CALL section_vals_val_get(rtp_section, "EPS_ITER", &
    3451          324 :                                 r_val=dft_control%rtp_control%eps_ener)
    3452              :       CALL section_vals_val_get(rtp_section, "INITIAL_WFN", &
    3453          324 :                                 i_val=dft_control%rtp_control%initial_wfn)
    3454              :       CALL section_vals_val_get(rtp_section, "HFX_BALANCE_IN_CORE", &
    3455          324 :                                 l_val=dft_control%rtp_control%hfx_redistribute)
    3456              :       CALL section_vals_val_get(rtp_section, "APPLY_WFN_MIX_INIT_RESTART", &
    3457          324 :                                 l_val=dft_control%rtp_control%apply_wfn_mix_init_restart)
    3458              :       CALL section_vals_val_get(rtp_section, "APPLY_DELTA_PULSE", &
    3459          324 :                                 l_val=dft_control%rtp_control%apply_delta_pulse)
    3460              :       CALL section_vals_val_get(rtp_section, "APPLY_DELTA_PULSE_MAG", &
    3461          324 :                                 l_val=dft_control%rtp_control%apply_delta_pulse_mag)
    3462              :       CALL section_vals_val_get(rtp_section, "VELOCITY_GAUGE", &
    3463          324 :                                 l_val=dft_control%rtp_control%velocity_gauge)
    3464              :       CALL section_vals_val_get(rtp_section, "VG_COM_NL", &
    3465          324 :                                 l_val=dft_control%rtp_control%nl_gauge_transform)
    3466              :       CALL section_vals_val_get(rtp_section, "PERIODIC", &
    3467          324 :                                 l_val=dft_control%rtp_control%periodic)
    3468              :       CALL section_vals_val_get(rtp_section, "DENSITY_PROPAGATION", &
    3469          324 :                                 l_val=dft_control%rtp_control%linear_scaling)
    3470              :       CALL section_vals_val_get(rtp_section, "MCWEENY_MAX_ITER", &
    3471          324 :                                 i_val=dft_control%rtp_control%mcweeny_max_iter)
    3472              :       CALL section_vals_val_get(rtp_section, "ACCURACY_REFINEMENT", &
    3473          324 :                                 i_val=dft_control%rtp_control%acc_ref)
    3474              :       CALL section_vals_val_get(rtp_section, "MCWEENY_EPS", &
    3475          324 :                                 r_val=dft_control%rtp_control%mcweeny_eps)
    3476              :       CALL section_vals_val_get(rtp_section, "DELTA_PULSE_SCALE", &
    3477          324 :                                 r_val=dft_control%rtp_control%delta_pulse_scale)
    3478              :       CALL section_vals_val_get(rtp_section, "DELTA_PULSE_DIRECTION", &
    3479          324 :                                 i_vals=tmp)
    3480         1296 :       dft_control%rtp_control%delta_pulse_direction = tmp
    3481              :       CALL section_vals_val_get(rtp_section, "SC_CHECK_START", &
    3482          324 :                                 i_val=dft_control%rtp_control%sc_check_start)
    3483          324 :       proj_mo_section => section_vals_get_subs_vals(rtp_section, "PRINT%PROJECTION_MO")
    3484          324 :       CALL section_vals_get(proj_mo_section, explicit=is_present)
    3485          324 :       IF (is_present) THEN
    3486            4 :          IF (dft_control%rtp_control%linear_scaling) THEN
    3487              :             CALL cp_abort(__LOCATION__, &
    3488              :                           "You have defined a time dependent projection of mos, but "// &
    3489              :                           "only the density matrix is propagated (DENSITY_PROPAGATION "// &
    3490              :                           ".TRUE.). Please either use MO-based real time DFT or do not "// &
    3491            0 :                           "define any PRINT%PROJECTION_MO section")
    3492              :          END IF
    3493            4 :          dft_control%rtp_control%is_proj_mo = .TRUE.
    3494              :       ELSE
    3495          320 :          dft_control%rtp_control%is_proj_mo = .FALSE.
    3496              :       END IF
    3497              :       ! Moment trace
    3498              :       local_moment_possible = (dft_control%rtp_control%rtp_method == rtp_method_bse .OR. &
    3499              :                                dft_control%rtp_control%rtp_method == rtp_method_bse_linearized) .OR. &
    3500          324 :                               ((.NOT. dft_control%rtp_control%periodic) .AND. dft_control%rtp_control%linear_scaling)
    3501              :       ! TODO : Implement for other moment operators
    3502          324 :       subsection => section_vals_get_subs_vals(rtp_section, "PRINT%MOMENTS")
    3503          324 :       CALL section_vals_get(subsection, explicit=is_present)
    3504              :       ! Trigger the flag
    3505              :       dft_control%rtp_control%save_local_moments = &
    3506          324 :          is_present .OR. dft_control%rtp_control%save_local_moments
    3507          324 :       IF (is_present .AND. (.NOT. local_moment_possible)) THEN
    3508              :          CALL cp_abort(__LOCATION__, "Moments trace printing only "// &
    3509              :                        "implemented in non-periodic systems in linear scaling. "// &
    3510            0 :                        "Please use DFT%PRINT%MOMENTS for other printing.")
    3511              :       END IF
    3512              :       CALL section_vals_val_get(rtp_section, "PRINT%MOMENTS%REFERENCE", &
    3513          324 :                                 i_val=dft_control%rtp_control%moment_trace_ref_type)
    3514              :       CALL section_vals_val_get(rtp_section, "PRINT%MOMENTS%REFERENCE_POINT", &
    3515          324 :                                 r_vals=dft_control%rtp_control%moment_trace_user_ref_point)
    3516              :       ! Moment Fourier transform
    3517          324 :       subsection => section_vals_get_subs_vals(rtp_section, "PRINT%MOMENTS_FT")
    3518          324 :       CALL section_vals_get(subsection, explicit=is_present)
    3519              :       ! Trigger the flag
    3520              :       dft_control%rtp_control%save_local_moments = &
    3521          324 :          is_present .OR. dft_control%rtp_control%save_local_moments
    3522          324 :       IF (is_present .AND. (.NOT. local_moment_possible)) THEN
    3523              :          ! Not implemented
    3524              :          CALL cp_abort(__LOCATION__, "Moments Fourier transform printing "// &
    3525            0 :                        "implemented only for non-periodic systems in linear scaling.")
    3526              :       END IF
    3527              :       ! General FT settings
    3528              :       CALL section_vals_val_get(rtp_section, "FT%DAMPING", &
    3529          324 :                                 r_val=dft_control%rtp_control%ft_damping)
    3530              :       CALL section_vals_val_get(rtp_section, "FT%START_TIME", &
    3531          324 :                                 r_val=dft_control%rtp_control%ft_t0)
    3532              :       ! Padé settings
    3533          324 :       subsection => section_vals_get_subs_vals(rtp_section, "FT%PADE")
    3534              :       CALL section_vals_val_get(subsection, "_SECTION_PARAMETERS_", &
    3535          324 :                                 l_val=dft_control%rtp_control%pade_requested)
    3536              :       CALL section_vals_val_get(subsection, "E_MIN", &
    3537          324 :                                 r_val=dft_control%rtp_control%pade_e_min)
    3538              :       CALL section_vals_val_get(subsection, "E_STEP", &
    3539          324 :                                 r_val=dft_control%rtp_control%pade_e_step)
    3540              :       CALL section_vals_val_get(subsection, "E_MAX", &
    3541          324 :                                 r_val=dft_control%rtp_control%pade_e_max)
    3542              :       CALL section_vals_val_get(subsection, "FIT_E_MIN", &
    3543          324 :                                 r_val=dft_control%rtp_control%pade_fit_e_min)
    3544              :       CALL section_vals_val_get(subsection, "FIT_E_MAX", &
    3545          324 :                                 r_val=dft_control%rtp_control%pade_fit_e_max)
    3546              :       ! If default settings used for fit_e_min/max, rewrite with appropriate values
    3547          324 :       IF (dft_control%rtp_control%pade_fit_e_min < 0) THEN
    3548          324 :          dft_control%rtp_control%pade_fit_e_min = dft_control%rtp_control%pade_e_min
    3549              :       END IF
    3550          324 :       IF (dft_control%rtp_control%pade_fit_e_max < 0) THEN
    3551          324 :          dft_control%rtp_control%pade_fit_e_max = dft_control%rtp_control%pade_e_max
    3552              :       END IF
    3553              :       ! Polarizability settings
    3554          324 :       subsection => section_vals_get_subs_vals(rtp_section, "PRINT%POLARIZABILITY")
    3555          324 :       CALL section_vals_get(subsection, explicit=is_present)
    3556              :       ! Trigger the flag
    3557              :       dft_control%rtp_control%save_local_moments = &
    3558          324 :          is_present .OR. dft_control%rtp_control%save_local_moments
    3559          324 :       IF (is_present .AND. (.NOT. local_moment_possible)) THEN
    3560              :          ! Not implemented
    3561              :          CALL cp_abort(__LOCATION__, "Polarizability printing "// &
    3562            0 :                        "implemented only for non-periodic systems.")
    3563              :       END IF
    3564          324 :       CALL section_vals_val_get(subsection, "ELEMENT", explicit=is_present, n_rep_val=n_elems)
    3565          324 :       NULLIFY (dft_control%rtp_control%print_pol_elements)
    3566          324 :       IF (is_present) THEN
    3567              :          ! Explicit list of elements
    3568              :          ! Allocate the array
    3569            0 :          ALLOCATE (dft_control%rtp_control%print_pol_elements(n_elems, 2))
    3570            0 :          DO i = 1, n_elems
    3571            0 :             CALL section_vals_val_get(subsection, "ELEMENT", i_vals=tmp, i_rep_val=i)
    3572            0 :             dft_control%rtp_control%print_pol_elements(i, :) = tmp(:)
    3573              :          END DO
    3574              :          ! Do basic sanity checks for pol_element
    3575            0 :          DO i = 1, n_elems
    3576            0 :             DO j = 1, 2
    3577            0 :                IF (dft_control%rtp_control%print_pol_elements(i, j) > 3 .OR. &
    3578            0 :                    dft_control%rtp_control%print_pol_elements(i, j) < 1) THEN
    3579            0 :                   CPABORT("Polarisation tensor element not 1,2 or 3 in at least one index")
    3580              :                END IF
    3581              :             END DO
    3582              :          END DO
    3583              :       END IF
    3584              : 
    3585              :       ! Finally, allow printing of FT observables also in the case when they are not explicitly
    3586              :       ! required, but they are available, i.e. non-periodic linear scaling calculation
    3587              :       dft_control%rtp_control%save_local_moments = &
    3588              :          dft_control%rtp_control%save_local_moments .OR. &
    3589          324 :          ((.NOT. dft_control%rtp_control%periodic) .AND. dft_control%rtp_control%linear_scaling)
    3590              : 
    3591          324 :    END SUBROUTINE read_rtp_section
    3592              : ! **************************************************************************************************
    3593              : !> \brief Tries to guess the elements of polarization to print
    3594              : !> \param dftc DFT parameters
    3595              : !> \param elems 2D array, where the guessed element indeces are stored
    3596              : !> \date 11.2025
    3597              : !> \author Stepan Marek
    3598              : ! **************************************************************************************************
    3599           90 :    SUBROUTINE guess_pol_elements(dftc, elems)
    3600              :       TYPE(dft_control_type)                             :: dftc
    3601              :       INTEGER, DIMENSION(:, :), POINTER                  :: elems
    3602              : 
    3603              :       INTEGER                                            :: i, i_nonzero, n_nonzero
    3604              :       LOGICAL                                            :: pol_vector_known
    3605              :       REAL(kind=dp), DIMENSION(3)                        :: pol_vector
    3606              : 
    3607           90 :       pol_vector_known = .FALSE.
    3608              : 
    3609              :       ! TODO : More relevant elements for magnetic pulse?
    3610           90 :       IF (dftc%rtp_control%apply_delta_pulse .OR. dftc%rtp_control%apply_delta_pulse_mag) THEN
    3611          336 :          pol_vector(:) = REAL(dftc%rtp_control%delta_pulse_direction(:), kind=dp)
    3612              :       ELSE
    3613              :          ! Maybe RT field is applied?
    3614           24 :          pol_vector(:) = dftc%efield_fields(1)%efield%polarisation(:)
    3615              :       END IF
    3616          360 :       IF (DOT_PRODUCT(pol_vector, pol_vector) > 0.0_dp) pol_vector_known = .TRUE.
    3617              : 
    3618              :       IF (.NOT. pol_vector_known) THEN
    3619            0 :          CPABORT("Cannot guess polarization elements - please specify!")
    3620              :       ELSE
    3621              :          ! Check whether just one element is non-zero
    3622              :          n_nonzero = 0
    3623          360 :          DO i = 1, 3
    3624          360 :             IF (pol_vector(i) /= 0.0_dp) THEN
    3625           90 :                n_nonzero = n_nonzero + 1
    3626           90 :                i_nonzero = i
    3627              :             END IF
    3628              :          END DO
    3629           90 :          IF (n_nonzero > 1) THEN
    3630              :             CALL cp_abort(__LOCATION__, &
    3631              :                           "More than one non-zero field elements - "// &
    3632            0 :                           "cannot guess polarizability elements - please specify!")
    3633           90 :          ELSE IF (n_nonzero == 0) THEN
    3634              :             CALL cp_abort(__LOCATION__, &
    3635              :                           "No non-zero field elements - "// &
    3636            0 :                           "cannot guess polarizability elements - please specify!")
    3637              :          ELSE
    3638              :             ! Clear guess can be made
    3639              :             NULLIFY (elems)
    3640           90 :             ALLOCATE (elems(3, 2))
    3641          360 :             DO i = 1, 3
    3642          270 :                elems(i, 1) = i
    3643          360 :                elems(i, 2) = i_nonzero
    3644              :             END DO
    3645              :          END IF
    3646              :       END IF
    3647           90 :    END SUBROUTINE guess_pol_elements
    3648              : 
    3649              : ! **************************************************************************************************
    3650              : !> \brief Parses the BLOCK_LIST keywords from the ADMM section
    3651              : !> \param admm_control ...
    3652              : !> \param dft_section ...
    3653              : ! **************************************************************************************************
    3654          524 :    SUBROUTINE read_admm_block_list(admm_control, dft_section)
    3655              :       TYPE(admm_control_type), POINTER                   :: admm_control
    3656              :       TYPE(section_vals_type), POINTER                   :: dft_section
    3657              : 
    3658              :       INTEGER                                            :: irep, list_size, n_rep
    3659          524 :       INTEGER, DIMENSION(:), POINTER                     :: tmplist
    3660              : 
    3661          524 :       NULLIFY (tmplist)
    3662              : 
    3663              :       CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%BLOCK_LIST", &
    3664          524 :                                 n_rep_val=n_rep)
    3665              : 
    3666         1102 :       ALLOCATE (admm_control%blocks(n_rep))
    3667              : 
    3668          560 :       DO irep = 1, n_rep
    3669              :          CALL section_vals_val_get(dft_section, "AUXILIARY_DENSITY_MATRIX_METHOD%BLOCK_LIST", &
    3670           36 :                                    i_rep_val=irep, i_vals=tmplist)
    3671           36 :          list_size = SIZE(tmplist)
    3672          108 :          ALLOCATE (admm_control%blocks(irep)%list(list_size))
    3673          732 :          admm_control%blocks(irep)%list(:) = tmplist(:)
    3674              :       END DO
    3675              : 
    3676          524 :    END SUBROUTINE read_admm_block_list
    3677              : 
    3678              : ! **************************************************************************************************
    3679              : !> \brief ...
    3680              : !> \param dft_control ...
    3681              : !> \param hairy_probes_section ...
    3682              : !> \param
    3683              : !> \param
    3684              : ! **************************************************************************************************
    3685            4 :    SUBROUTINE read_hairy_probes_sections(dft_control, hairy_probes_section)
    3686              :       TYPE(dft_control_type), POINTER                    :: dft_control
    3687              :       TYPE(section_vals_type), POINTER                   :: hairy_probes_section
    3688              : 
    3689              :       INTEGER                                            :: i, j, jj, kk, n_rep
    3690            4 :       INTEGER, DIMENSION(:), POINTER                     :: tmplist
    3691              : 
    3692           12 :       DO i = 1, SIZE(dft_control%probe)
    3693            8 :          NULLIFY (dft_control%probe(i)%atom_ids)
    3694              : 
    3695            8 :          CALL section_vals_val_get(hairy_probes_section, "ATOM_IDS", i_rep_section=i, n_rep_val=n_rep)
    3696            8 :          jj = 0
    3697           16 :          DO kk = 1, n_rep
    3698            8 :             CALL section_vals_val_get(hairy_probes_section, "ATOM_IDS", i_rep_section=i, i_rep_val=kk, i_vals=tmplist)
    3699           16 :             jj = jj + SIZE(tmplist)
    3700              :          END DO
    3701              : 
    3702            8 :          dft_control%probe(i)%natoms = jj
    3703            8 :          IF (dft_control%probe(i)%natoms < 1) THEN
    3704            0 :             CPABORT("Need at least 1 atom to use hair probes formalism")
    3705              :          END IF
    3706           24 :          ALLOCATE (dft_control%probe(i)%atom_ids(dft_control%probe(i)%natoms))
    3707              : 
    3708            8 :          jj = 0
    3709           16 :          DO kk = 1, n_rep
    3710            8 :             CALL section_vals_val_get(hairy_probes_section, "ATOM_IDS", i_rep_section=i, i_rep_val=kk, i_vals=tmplist)
    3711           24 :             DO j = 1, SIZE(tmplist)
    3712            8 :                jj = jj + 1
    3713           16 :                dft_control%probe(i)%atom_ids(jj) = tmplist(j)
    3714              :             END DO
    3715              :          END DO
    3716              : 
    3717            8 :          CALL section_vals_val_get(hairy_probes_section, "MU", i_rep_section=i, r_val=dft_control%probe(i)%mu)
    3718              : 
    3719            8 :          CALL section_vals_val_get(hairy_probes_section, "T", i_rep_section=i, r_val=dft_control%probe(i)%T)
    3720              : 
    3721            8 :          CALL section_vals_val_get(hairy_probes_section, "ALPHA", i_rep_section=i, r_val=dft_control%probe(i)%alpha)
    3722              : 
    3723           20 :          CALL section_vals_val_get(hairy_probes_section, "eps_hp", i_rep_section=i, r_val=dft_control%probe(i)%eps_hp)
    3724              :       END DO
    3725              : 
    3726            4 :    END SUBROUTINE read_hairy_probes_sections
    3727              : ! **************************************************************************************************
    3728              : 
    3729              : END MODULE cp_control_utils
        

Generated by: LCOV version 2.0-1