LCOV - code coverage report
Current view: top level - src - cp_control_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:06f838d) Lines: 87.3 % 1838 1604
Test Date: 2026-06-05 07:04:50 Functions: 100.0 % 17 17

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

Generated by: LCOV version 2.0-1