LCOV - code coverage report
Current view: top level - src - cp_control_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:936074a) Lines: 88.0 % 1594 1402
Test Date: 2025-12-04 06:27:48 Functions: 100.0 % 15 15

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

Generated by: LCOV version 2.0-1