!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2020 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!

! **************************************************************************************************
!> \par History
!>      - Merged with the Quickstep MODULE method_specification (17.01.2002,MK)
!>      - USE statements cleaned, added
!>        (25.09.2002,MK)
!>      - Added more LSD structure (01.2003,Joost VandeVondele)
!>      - New molecule data types introduced (Sep. 2003,MK)
!>      - Cleaning; getting rid of pnode (02.10.2003,MK)
!>      - Sub-system setup added (08.10.2003,MK)
!> \author MK (18.05.2000)
! **************************************************************************************************
MODULE qs_environment
   USE almo_scf_env_methods,            ONLY: almo_scf_env_create
   USE atom_kind_orbitals,              ONLY: calculate_atomic_relkin
   USE atomic_kind_types,               ONLY: atomic_kind_type
   USE auto_basis,                      ONLY: create_aux_fit_basis_set,&
                                              create_lri_aux_basis_set,&
                                              create_ri_aux_basis_set
   USE basis_set_container_types,       ONLY: add_basis_set_to_container
   USE basis_set_types,                 ONLY: basis_sort_zet,&
                                              get_gto_basis_set,&
                                              gto_basis_set_type
   USE bibliography,                    ONLY: Iannuzzi2006,&
                                              Iannuzzi2007,&
                                              cite_reference,&
                                              cp2kqs2020
   USE cell_types,                      ONLY: cell_type
   USE cp_blacs_env,                    ONLY: cp_blacs_env_create,&
                                              cp_blacs_env_release,&
                                              cp_blacs_env_retain,&
                                              cp_blacs_env_type
   USE cp_control_types,                ONLY: dft_control_release,&
                                              dft_control_type,&
                                              dftb_control_type,&
                                              gapw_control_type,&
                                              qs_control_type,&
                                              semi_empirical_control_type,&
                                              xtb_control_type
   USE cp_control_utils,                ONLY: read_ddapc_section,&
                                              read_dft_control,&
                                              read_mgrid_section,&
                                              read_qs_section,&
                                              read_tddfpt2_control,&
                                              read_tddfpt_control,&
                                              write_dft_control,&
                                              write_qs_control
   USE cp_ddapc_types,                  ONLY: cp_ddapc_ewald_create
   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
                                              cp_logger_get_default_io_unit,&
                                              cp_logger_type
   USE cp_output_handling,              ONLY: cp_p_file,&
                                              cp_print_key_finished_output,&
                                              cp_print_key_should_output,&
                                              cp_print_key_unit_nr
   USE cp_para_env,                     ONLY: cp_para_env_retain
   USE cp_para_types,                   ONLY: cp_para_env_type
   USE cp_subsys_types,                 ONLY: cp_subsys_type
   USE cp_symmetry,                     ONLY: write_symmetry
   USE distribution_1d_types,           ONLY: distribution_1d_release,&
                                              distribution_1d_type
   USE distribution_methods,            ONLY: distribute_molecules_1d
   USE dm_ls_scf_create,                ONLY: ls_scf_create
   USE ec_env_types,                    ONLY: energy_correction_type
   USE ec_environment,                  ONLY: ec_env_create
   USE et_coupling_types,               ONLY: et_coupling_create
   USE ewald_environment_types,         ONLY: ewald_env_create,&
                                              ewald_env_get,&
                                              ewald_env_release,&
                                              ewald_env_set,&
                                              ewald_environment_type,&
                                              read_ewald_section,&
                                              read_ewald_section_tb
   USE ewald_pw_methods,                ONLY: ewald_pw_grid_update
   USE ewald_pw_types,                  ONLY: ewald_pw_create,&
                                              ewald_pw_release,&
                                              ewald_pw_type
   USE external_potential_types,        ONLY: get_potential,&
                                              init_potential,&
                                              set_potential
   USE fist_nonbond_env_types,          ONLY: fist_nonbond_env_create,&
                                              fist_nonbond_env_type
   USE gamma,                           ONLY: init_md_ftable
   USE global_types,                    ONLY: global_environment_type
   USE hartree_local_methods,           ONLY: init_coulomb_local
   USE header,                          ONLY: dftb_header,&
                                              qs_header,&
                                              se_header,&
                                              xtb_header
   USE hfx_types,                       ONLY: hfx_create
   USE input_constants,                 ONLY: &
        dispersion_d3, do_et_ddapc, do_method_am1, do_method_dftb, do_method_gapw, &
        do_method_gapw_xc, do_method_gpw, do_method_lrigpw, do_method_mndo, do_method_mndod, &
        do_method_ofgpw, do_method_pdg, do_method_pm3, do_method_pm6, do_method_pm6fm, &
        do_method_pnnl, do_method_rigpw, do_method_rm1, do_method_xtb, do_qmmm_gauss, &
        do_qmmm_swave, general_roks, kg_tnadd_embed_ri, rel_none, rel_trans_atom, &
        vdw_pairpot_dftd3, vdw_pairpot_dftd3bj, wfi_aspc_nr, wfi_linear_ps_method_nr, &
        wfi_linear_wf_method_nr, wfi_ps_method_nr, wfi_use_guess_method_nr, xc_vdw_fun_none, &
        xc_vdw_fun_nonloc, xc_vdw_fun_pairpot
   USE input_section_types,             ONLY: section_vals_get,&
                                              section_vals_get_subs_vals,&
                                              section_vals_type,&
                                              section_vals_val_get
   USE kg_environment,                  ONLY: kg_env_create
   USE kinds,                           ONLY: dp
   USE kpoint_methods,                  ONLY: kpoint_env_initialize,&
                                              kpoint_initialize,&
                                              kpoint_initialize_mos
   USE kpoint_types,                    ONLY: kpoint_create,&
                                              kpoint_type,&
                                              read_kpoint_section,&
                                              write_kpoint_info
   USE lri_environment_init,            ONLY: lri_env_basis,&
                                              lri_env_init
   USE lri_environment_types,           ONLY: lri_environment_type
   USE machine,                         ONLY: m_flush
   USE mathconstants,                   ONLY: pi
   USE molecule_kind_types,             ONLY: molecule_kind_type,&
                                              write_molecule_kind_set
   USE molecule_types,                  ONLY: molecule_type
   USE mp2_setup,                       ONLY: read_mp2_section
   USE mp2_types,                       ONLY: mp2_env_create,&
                                              mp2_type
   USE multipole_types,                 ONLY: do_multipole_none
   USE orbital_pointers,                ONLY: init_orbital_pointers
   USE orbital_transformation_matrices, ONLY: init_spherical_harmonics
   USE particle_methods,                ONLY: write_particle_distances,&
                                              write_qs_particle_coordinates,&
                                              write_structure_data
   USE particle_types,                  ONLY: particle_type
   USE pw_env_types,                    ONLY: pw_env_type
   USE qmmm_types_low,                  ONLY: qmmm_env_qm_type
   USE qs_dftb_parameters,              ONLY: qs_dftb_param_init
   USE qs_dftb_types,                   ONLY: qs_dftb_pairpot_type
   USE qs_dispersion_nonloc,            ONLY: qs_dispersion_nonloc_init
   USE qs_dispersion_pairpot,           ONLY: qs_dispersion_pairpot_init
   USE qs_dispersion_types,             ONLY: qs_dispersion_type
   USE qs_dispersion_utils,             ONLY: qs_dispersion_env_set,&
                                              qs_write_dispersion
   USE qs_energy_types,                 ONLY: allocate_qs_energy,&
                                              qs_energy_type
   USE qs_environment_methods,          ONLY: qs_env_setup
   USE qs_environment_types,            ONLY: get_qs_env,&
                                              qs_environment_type,&
                                              set_qs_env
   USE qs_force_types,                  ONLY: qs_force_type
   USE qs_gcp_types,                    ONLY: qs_gcp_type
   USE qs_gcp_utils,                    ONLY: qs_gcp_env_set,&
                                              qs_gcp_init
   USE qs_interactions,                 ONLY: init_interaction_radii,&
                                              init_se_nlradius,&
                                              write_core_charge_radii,&
                                              write_paw_radii,&
                                              write_pgf_orb_radii,&
                                              write_ppl_radii,&
                                              write_ppnl_radii
   USE qs_kind_types,                   ONLY: &
        check_qs_kind_set, get_qs_kind, get_qs_kind_set, init_gapw_basis_set, init_gapw_nlcc, &
        init_qs_kind_set, qs_kind_type, set_qs_kind, write_gto_basis_sets, write_qs_kind_set
   USE qs_ks_types,                     ONLY: qs_ks_env_create,&
                                              qs_ks_env_type,&
                                              qs_ks_release,&
                                              set_ks_env
   USE qs_local_rho_types,              ONLY: local_rho_type
   USE qs_mo_types,                     ONLY: allocate_mo_set,&
                                              mo_set_p_type
   USE qs_rho0_ggrid,                   ONLY: rho0_s_grid_create
   USE qs_rho0_methods,                 ONLY: init_rho0
   USE qs_rho0_types,                   ONLY: rho0_mpole_type
   USE qs_rho_atom_methods,             ONLY: init_rho_atom
   USE qs_rho_atom_types,               ONLY: rho_atom_type
   USE qs_subsys_methods,               ONLY: qs_subsys_create
   USE qs_subsys_types,                 ONLY: qs_subsys_get,&
                                              qs_subsys_release,&
                                              qs_subsys_set,&
                                              qs_subsys_type
   USE qs_wf_history_methods,           ONLY: wfi_create,&
                                              wfi_create_for_kp
   USE qs_wf_history_types,             ONLY: qs_wf_history_type,&
                                              wfi_release
   USE rel_control_types,               ONLY: rel_c_create,&
                                              rel_c_read_parameters,&
                                              rel_c_release,&
                                              rel_control_type
   USE scf_control_types,               ONLY: scf_c_create,&
                                              scf_c_read_parameters,&
                                              scf_c_release,&
                                              scf_c_write_parameters,&
                                              scf_control_type
   USE semi_empirical_expns3_methods,   ONLY: semi_empirical_expns3_setup
   USE semi_empirical_int_arrays,       ONLY: init_se_intd_array
   USE semi_empirical_mpole_methods,    ONLY: nddo_mpole_setup
   USE semi_empirical_mpole_types,      ONLY: nddo_mpole_type
   USE semi_empirical_store_int_types,  ONLY: semi_empirical_si_create,&
                                              semi_empirical_si_type
   USE semi_empirical_types,            ONLY: se_taper_create,&
                                              se_taper_type
   USE semi_empirical_utils,            ONLY: se_cutoff_compatible
   USE transport,                       ONLY: transport_env_create
   USE xtb_parameters,                  ONLY: init_xtb_basis,&
                                              xtb_parameters_init,&
                                              xtb_parameters_read,&
                                              xtb_parameters_set
   USE xtb_types,                       ONLY: allocate_xtb_atom_param
#include "./base/base_uses.f90"

   IMPLICIT NONE

   PRIVATE

   ! *** Global parameters ***
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_environment'

   ! *** Public subroutines ***
   PUBLIC :: qs_init

CONTAINS

! **************************************************************************************************
!> \brief Read the input and the database files for the setup of the
!>      QUICKSTEP environment.
!> \param qs_env ...
!> \param para_env ...
!> \param root_section ...
!> \param globenv ...
!> \param cp_subsys ...
!> \param kpoint_env ...
!> \param cell ...
!> \param cell_ref ...
!> \param qmmm ...
!> \param qmmm_env_qm ...
!> \param force_env_section ...
!> \param subsys_section ...
!> \param use_motion_section ...
!> \author Creation (22.05.2000,MK)
! **************************************************************************************************
   SUBROUTINE qs_init(qs_env, para_env, root_section, globenv, cp_subsys, kpoint_env, cell, cell_ref, &
                      qmmm, qmmm_env_qm, force_env_section, subsys_section, use_motion_section)

      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(cp_para_env_type), POINTER                    :: para_env
      TYPE(section_vals_type), OPTIONAL, POINTER         :: root_section
      TYPE(global_environment_type), OPTIONAL, POINTER   :: globenv
      TYPE(cp_subsys_type), OPTIONAL, POINTER            :: cp_subsys
      TYPE(kpoint_type), OPTIONAL, POINTER               :: kpoint_env
      TYPE(cell_type), OPTIONAL, POINTER                 :: cell, cell_ref
      LOGICAL, INTENT(IN), OPTIONAL                      :: qmmm
      TYPE(qmmm_env_qm_type), OPTIONAL, POINTER          :: qmmm_env_qm
      TYPE(section_vals_type), POINTER                   :: force_env_section, subsys_section
      LOGICAL, INTENT(IN)                                :: use_motion_section

      INTEGER                                            :: ikind, method_id, nelectron_total, nkind
      LOGICAL                                            :: do_admm_rpa, do_et, do_exx, do_hfx, &
                                                            do_kpoints, is_semi, mp2_present, &
                                                            my_qmmm, qmmm_decoupl, silent, &
                                                            use_ref_cell
      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: rtmat
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cell_type), POINTER                           :: my_cell, my_cell_ref
      TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(energy_correction_type), POINTER              :: ec_env
      TYPE(kpoint_type), POINTER                         :: kpoints
      TYPE(lri_environment_type), POINTER                :: lri_env
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      TYPE(qs_ks_env_type), POINTER                      :: ks_env
      TYPE(qs_subsys_type), POINTER                      :: subsys
      TYPE(qs_wf_history_type), POINTER                  :: wf_history
      TYPE(rel_control_type), POINTER                    :: rel_control
      TYPE(scf_control_type), POINTER                    :: scf_control
      TYPE(section_vals_type), POINTER                   :: dft_section, et_coupling_section, &
                                                            hfx_section, kpoint_section, &
                                                            mp2_section, transport_section

      NULLIFY (my_cell, my_cell_ref, atomic_kind_set, particle_set, &
               qs_kind_set, kpoint_section, dft_section, &
               subsys, ks_env, dft_control, blacs_env)

      CALL set_qs_env(qs_env, input=force_env_section)
      IF (.NOT. ASSOCIATED(subsys_section)) THEN
         subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
      END IF

      ! QMMM
      my_qmmm = .FALSE.
      IF (PRESENT(qmmm)) my_qmmm = qmmm
      qmmm_decoupl = .FALSE.
      IF (PRESENT(qmmm_env_qm)) THEN
         IF (qmmm_env_qm%qmmm_coupl_type == do_qmmm_gauss .OR. &
             qmmm_env_qm%qmmm_coupl_type == do_qmmm_swave) THEN
            ! For GAUSS/SWAVE methods there could be a DDAPC decoupling requested
            qmmm_decoupl = my_qmmm .AND. qmmm_env_qm%periodic .AND. qmmm_env_qm%multipole
         END IF
         qs_env%qmmm_env_qm => qmmm_env_qm
      END IF
      CALL set_qs_env(qs_env=qs_env, qmmm=my_qmmm)

      ! Possibly initialize arrays for SE
      CALL section_vals_val_get(force_env_section, "DFT%QS%METHOD", i_val=method_id)
      SELECT CASE (method_id)
      CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pdg, &
            do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
         CALL init_se_intd_array()
         is_semi = .TRUE.
      CASE (do_method_xtb, do_method_dftb)
         is_semi = .TRUE.
      CASE DEFAULT
         is_semi = .FALSE.
      END SELECT

      CALL qs_subsys_create(subsys, para_env, &
                            force_env_section=force_env_section, &
                            subsys_section=subsys_section, &
                            use_motion_section=use_motion_section, &
                            root_section=root_section, &
                            cp_subsys=cp_subsys, cell=cell, cell_ref=cell_ref, &
                            elkind=is_semi)

      CALL qs_ks_env_create(ks_env)
      CALL set_ks_env(ks_env, subsys=subsys)
      CALL set_qs_env(qs_env, ks_env=ks_env)

      CALL qs_subsys_get(subsys, &
                         cell=my_cell, &
                         cell_ref=my_cell_ref, &
                         use_ref_cell=use_ref_cell, &
                         atomic_kind_set=atomic_kind_set, &
                         qs_kind_set=qs_kind_set, &
                         particle_set=particle_set)

      CALL set_ks_env(ks_env, para_env=para_env)
      IF (PRESENT(globenv)) THEN
         CALL cp_blacs_env_create(blacs_env, para_env, globenv%blacs_grid_layout, &
                                  globenv%blacs_repeatable)
      ELSE
         CALL cp_blacs_env_create(blacs_env, para_env)
      END IF
      CALL set_ks_env(ks_env, blacs_env=blacs_env)
      CALL cp_blacs_env_release(blacs_env)

      !   *** Setup the grids for the G-space Interpolation if any
      CALL cp_ddapc_ewald_create(qs_env%cp_ddapc_ewald, qmmm_decoupl, my_cell, &
                                 force_env_section, subsys_section, para_env)

      ! kpoints
      IF (PRESENT(kpoint_env)) THEN
         silent = .TRUE.
         kpoints => kpoint_env
         CALL set_qs_env(qs_env=qs_env, kpoints=kpoints)
         CALL kpoint_initialize(kpoints, particle_set, my_cell)
      ELSE
         silent = .FALSE.
         NULLIFY (kpoints)
         CALL kpoint_create(kpoints)
         CALL set_qs_env(qs_env=qs_env, kpoints=kpoints)
         kpoint_section => section_vals_get_subs_vals(qs_env%input, "DFT%KPOINTS")
         CALL read_kpoint_section(kpoints, kpoint_section, my_cell%hmat)
         CALL kpoint_initialize(kpoints, particle_set, my_cell)
         dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
         CALL write_kpoint_info(kpoints, dft_section)
      END IF

      CALL qs_init_subsys(qs_env, para_env, subsys, my_cell, my_cell_ref, use_ref_cell, &
                          subsys_section, silent=silent)

      CALL get_qs_env(qs_env, dft_control=dft_control)
      IF (method_id == do_method_lrigpw .OR. dft_control%qs_control%lri_optbas) THEN
         CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
         CALL lri_env_basis("LRI", qs_env, lri_env, qs_kind_set)
      ELSE IF (method_id == do_method_rigpw) THEN
         CALL cp_warn(__LOCATION__, "Experimental code: "// &
                      "RIGPW should only be used for testing.")
         CALL get_qs_env(qs_env=qs_env, lri_env=lri_env)
         CALL lri_env_basis("RI", qs_env, lri_env, qs_kind_set)
      END IF

      IF (my_qmmm .AND. PRESENT(qmmm_env_qm) .AND. .NOT. dft_control%qs_control%commensurate_mgrids) THEN
         IF (qmmm_env_qm%qmmm_coupl_type == do_qmmm_gauss .OR. qmmm_env_qm%qmmm_coupl_type == do_qmmm_swave) THEN
            CALL cp_abort(__LOCATION__, "QM/MM with coupling GAUSS or S-WAVE requires "// &
                          "keyword FORCE_EVAL/DFT/MGRID/COMMENSURATE to be enabled.")
         ENDIF
      ENDIF

      ! more kpoint stuff
      kpoints%para_env => para_env
      CALL cp_para_env_retain(para_env)
      CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env)
      kpoints%blacs_env_all => blacs_env
      CALL cp_blacs_env_retain(blacs_env)
      CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
      IF (do_kpoints) THEN
         CALL kpoint_env_initialize(kpoints)
         CALL kpoint_initialize_mos(kpoints, qs_env%mos)
         CALL get_qs_env(qs_env=qs_env, wf_history=wf_history)
         CALL wfi_create_for_kp(wf_history)
      END IF

      do_hfx = .FALSE.
      hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%HF")
      CALL section_vals_get(hfx_section, explicit=do_hfx)
      CALL get_qs_env(qs_env, dft_control=dft_control, scf_control=scf_control, nelectron_total=nelectron_total)
      IF (do_hfx) THEN
         ! Retrieve particle_set and atomic_kind_set (needed for both kinds of initialization)
         CALL hfx_create(qs_env%x_data, para_env, hfx_section, atomic_kind_set, &
                         qs_kind_set, particle_set, dft_control, my_cell, do_ot=scf_control%use_ot, &
                         nelectron_total=nelectron_total)
      END IF

      mp2_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION")
      CALL section_vals_get(mp2_section, explicit=mp2_present)
      IF (mp2_present) THEN
         CPASSERT(ASSOCIATED(qs_env%mp2_env))
         CALL read_mp2_section(qs_env%input, qs_env%mp2_env)
         ! create the EXX section if necessary
         do_exx = .FALSE.
         hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
         CALL section_vals_get(hfx_section, explicit=do_exx)
         IF (do_exx) THEN
            ! Retrieve particle_set and atomic_kind_set (needed for both kinds of initialization)

            ! do_exx in call of hfx_create decides whether to go without ADMM (do_exx=.TRUE.) or with
            ! ADMM (do_exx=.FALSE.)
            CALL section_vals_val_get(mp2_section, "RI_RPA%ADMM", l_val=do_admm_rpa)

            CALL hfx_create(qs_env%mp2_env%ri_rpa%x_data, para_env, hfx_section, atomic_kind_set, &
                            qs_kind_set, particle_set, dft_control, my_cell, do_exx=(.NOT. do_admm_rpa), &
                            do_ot=scf_control%use_ot, nelectron_total=nelectron_total)

         END IF
      END IF

      IF (dft_control%qs_control%do_kg) THEN
         CALL cite_reference(Iannuzzi2006)
         CALL kg_env_create(qs_env, qs_env%kg_env, qs_kind_set, qs_env%input)
      END IF

      dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
      CALL section_vals_val_get(dft_section, "ENERGY_CORRECTION%_SECTION_PARAMETERS_", &
                                l_val=qs_env%energy_correction)
      NULLIFY (ec_env)
      CALL ec_env_create(qs_env, ec_env, dft_section)
      CALL set_qs_env(qs_env, ec_env=ec_env)

      et_coupling_section => section_vals_get_subs_vals(qs_env%input, &
                                                        "PROPERTIES%ET_COUPLING")
      CALL section_vals_get(et_coupling_section, explicit=do_et)
      IF (do_et) CALL et_coupling_create(qs_env%et_coupling)

      transport_section => section_vals_get_subs_vals(qs_env%input, "DFT%TRANSPORT")
      CALL section_vals_get(transport_section, explicit=qs_env%do_transport)
      IF (qs_env%do_transport) THEN
         CALL transport_env_create(qs_env)
      END IF

      IF (dft_control%qs_control%do_ls_scf) THEN
         CALL ls_scf_create(qs_env)
      ENDIF

      IF (dft_control%qs_control%do_almo_scf) THEN
         CALL almo_scf_env_create(qs_env)
      ENDIF

      ! see if we have atomic relativistic corrections
      CALL get_qs_env(qs_env, rel_control=rel_control)
      IF (rel_control%rel_method /= rel_none) THEN
         IF (rel_control%rel_transformation == rel_trans_atom) THEN
            nkind = SIZE(atomic_kind_set)
            DO ikind = 1, nkind
               NULLIFY (rtmat)
               CALL calculate_atomic_relkin(atomic_kind_set(ikind), qs_kind_set(ikind), rel_control, rtmat)
               IF (ASSOCIATED(rtmat)) CALL set_qs_kind(qs_kind_set(ikind), reltmat=rtmat)
            END DO
         END IF
      END IF

      CALL qs_subsys_release(subsys)
      CALL qs_ks_release(ks_env)

   END SUBROUTINE qs_init

! **************************************************************************************************
!> \brief Initialize the qs environment (subsys)
!> \param qs_env ...
!> \param para_env ...
!> \param subsys ...
!> \param cell ...
!> \param cell_ref ...
!> \param use_ref_cell ...
!> \param subsys_section ...
!> \param silent ...
!> \author Creation (22.05.2000,MK)
! **************************************************************************************************
   SUBROUTINE qs_init_subsys(qs_env, para_env, subsys, cell, cell_ref, use_ref_cell, subsys_section, &
                             silent)

      TYPE(qs_environment_type), POINTER                 :: qs_env
      TYPE(cp_para_env_type), POINTER                    :: para_env
      TYPE(qs_subsys_type), POINTER                      :: subsys
      TYPE(cell_type), POINTER                           :: cell, cell_ref
      LOGICAL, INTENT(in)                                :: use_ref_cell
      TYPE(section_vals_type), POINTER                   :: subsys_section
      LOGICAL, INTENT(in), OPTIONAL                      :: silent

      CHARACTER(len=*), PARAMETER                        :: routineN = 'qs_init_subsys'

      CHARACTER(len=2)                                   :: element_symbol
      INTEGER :: handle, ikind, ispin, iw, lmax_sphere, maxl, maxlgto, maxlgto_lri, maxlppl, &
         maxlppnl, method_id, multiplicity, my_ival, n_ao, n_ao_aux_fit, n_mo_add, natom, &
         nelectron, ngauss, nkind, output_unit, sort_basis, tnadd_method
      INTEGER, DIMENSION(2)                              :: n_mo, nelectron_spin
      LOGICAL :: all_potential_present, be_silent, do_kpoints, do_ri_hfx, do_ri_mp2, do_ri_rpa, &
         do_ri_sos_mp2, do_wfc_im_time, e1terms, has_unit_metric, lribas, mp2_present, orb_gradient
      REAL(dp)                                           :: alpha, ccore, ewald_rcut, maxocc, &
                                                            verlet_skin, zeff_correction
      REAL(KIND=dp)                                      :: total_zeff_corr
      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      TYPE(cp_logger_type), POINTER                      :: logger
      TYPE(dft_control_type), POINTER                    :: dft_control
      TYPE(dftb_control_type), POINTER                   :: dftb_control
      TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
      TYPE(ewald_environment_type), POINTER              :: ewald_env
      TYPE(ewald_pw_type), POINTER                       :: ewald_pw
      TYPE(fist_nonbond_env_type), POINTER               :: se_nonbond_env
      TYPE(gapw_control_type), POINTER                   :: gapw_control
      TYPE(gto_basis_set_type), POINTER                  :: aux_fit_basis, lri_aux_basis, &
                                                            ri_aux_basis_set, ri_hfx_basis, &
                                                            ri_xas_basis, tmp_basis_set
      TYPE(local_rho_type), POINTER                      :: local_rho_set
      TYPE(lri_environment_type), POINTER                :: lri_env
      TYPE(mo_set_p_type), DIMENSION(:), POINTER         :: mos, mos_aux_fit, mos_last_converged
      TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
      TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
      TYPE(mp2_type), POINTER                            :: mp2_env
      TYPE(nddo_mpole_type), POINTER                     :: se_nddo_mpole
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(pw_env_type), POINTER                         :: pw_env
      TYPE(qs_control_type), POINTER                     :: qs_control
      TYPE(qs_dftb_pairpot_type), DIMENSION(:, :), &
         POINTER                                         :: dftb_potential
      TYPE(qs_dispersion_type), POINTER                  :: dispersion_env
      TYPE(qs_energy_type), POINTER                      :: energy
      TYPE(qs_force_type), DIMENSION(:), POINTER         :: force
      TYPE(qs_gcp_type), POINTER                         :: gcp_env
      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      TYPE(qs_kind_type), POINTER                        :: qs_kind
      TYPE(qs_ks_env_type), POINTER                      :: ks_env
      TYPE(qs_wf_history_type), POINTER                  :: wf_history
      TYPE(rel_control_type), POINTER                    :: rel_control
      TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole
      TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho_atom_set
      TYPE(scf_control_type), POINTER                    :: scf_control
      TYPE(se_taper_type), POINTER                       :: se_taper
      TYPE(section_vals_type), POINTER :: dft_section, et_coupling_section, et_ddapc_section, &
         ewald_section, lri_section, mp2_section, nl_section, poisson_section, pp_section, &
         print_section, qs_section, se_section, tddfpt_section, xc_section, xtb_section
      TYPE(semi_empirical_control_type), POINTER         :: se_control
      TYPE(semi_empirical_si_type), POINTER              :: se_store_int_env
      TYPE(xtb_control_type), POINTER                    :: xtb_control

      CALL timeset(routineN, handle)
      NULLIFY (logger)
      logger => cp_get_default_logger()
      output_unit = cp_logger_get_default_io_unit(logger)

      be_silent = .FALSE.
      IF (PRESENT(silent)) be_silent = silent

      CALL cite_reference(cp2kqs2020)

      ! Initialise the Quickstep environment
      NULLIFY (mos, se_taper, mos_aux_fit)
      NULLIFY (dft_control)
      NULLIFY (energy)
      NULLIFY (force)
      NULLIFY (local_molecules)
      NULLIFY (local_particles)
      NULLIFY (scf_control)
      NULLIFY (dft_section)
      NULLIFY (et_coupling_section)
      NULLIFY (ks_env)
      NULLIFY (mos_last_converged)
      dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
      qs_section => section_vals_get_subs_vals(dft_section, "QS")
      et_coupling_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%ET_COUPLING")
      ! reimplemented TDDFPT
      tddfpt_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%TDDFPT")

      CALL qs_subsys_get(subsys, particle_set=particle_set, &
                         qs_kind_set=qs_kind_set, &
                         atomic_kind_set=atomic_kind_set, &
                         molecule_set=molecule_set, &
                         molecule_kind_set=molecule_kind_set)

      !   *** Read the input section with the DFT control parameters ***
      CALL read_dft_control(dft_control, dft_section)

      ! original implementation of TDDFPT
      IF (dft_control%do_tddfpt_calculation) THEN
         CALL read_tddfpt_control(dft_control%tddfpt_control, dft_section)
      END IF

      !   *** Print the Quickstep program banner (copyright and version number) ***
      IF (.NOT. be_silent) THEN
         iw = cp_print_key_unit_nr(logger, dft_section, "PRINT%PROGRAM_BANNER", extension=".Log")
         CALL section_vals_val_get(qs_section, "METHOD", i_val=method_id)
         SELECT CASE (method_id)
         CASE DEFAULT
            CALL qs_header(iw)
         CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pdg, &
               do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
            CALL se_header(iw)
         CASE (do_method_dftb)
            CALL dftb_header(iw)
         CASE (do_method_xtb)
            CALL xtb_header(iw)
         END SELECT
         CALL cp_print_key_finished_output(iw, logger, dft_section, &
                                           "PRINT%PROGRAM_BANNER")
      END IF

      !   *** Read the input section with the Quickstep control parameters ***
      CALL read_qs_section(dft_control%qs_control, qs_section)
      CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
      IF (do_kpoints) THEN
         ! reset some of the settings for wfn extrapolation for kpoints
         SELECT CASE (dft_control%qs_control%wf_interpolation_method_nr)
         CASE (wfi_linear_wf_method_nr, wfi_linear_ps_method_nr, wfi_ps_method_nr, wfi_aspc_nr)
            dft_control%qs_control%wf_interpolation_method_nr = wfi_use_guess_method_nr
         END SELECT
      END IF

      !   *******  check if any kind of electron transfer calculation has to be performed
      CALL section_vals_val_get(et_coupling_section, "TYPE_OF_CONSTRAINT", i_val=my_ival)
      dft_control%qs_control%et_coupling_calc = .FALSE.
      IF (my_ival == do_et_ddapc) THEN
         et_ddapc_section => section_vals_get_subs_vals(et_coupling_section, "DDAPC_RESTRAINT_A")
         dft_control%qs_control%et_coupling_calc = .TRUE.
         dft_control%qs_control%ddapc_restraint = .TRUE.
         CALL read_ddapc_section(dft_control%qs_control, ddapc_restraint_section=et_ddapc_section)
      ENDIF

      CALL read_mgrid_section(dft_control%qs_control, dft_section)

      ! reimplemented TDDFPT
      CALL read_tddfpt2_control(dft_control%tddfpt2_control, tddfpt_section, dft_control%qs_control)

      !   Create relativistic control section
      CALL rel_c_create(rel_control)
      CALL rel_c_read_parameters(rel_control, dft_section)
      CALL set_qs_env(qs_env, rel_control=rel_control)
      CALL rel_c_release(rel_control)

      !   *** Read DFTB parameter files ***
      IF (dft_control%qs_control%method_id == do_method_dftb) THEN
         NULLIFY (ewald_env, ewald_pw, dftb_potential)
         dftb_control => dft_control%qs_control%dftb_control
         CALL qs_dftb_param_init(atomic_kind_set, qs_kind_set, dftb_control, dftb_potential, &
                                 subsys_section=subsys_section, para_env=para_env)
         CALL set_qs_env(qs_env, dftb_potential=dftb_potential)
         ! check for Ewald
         IF (dftb_control%do_ewald) THEN
            CALL ewald_env_create(ewald_env, para_env)
            poisson_section => section_vals_get_subs_vals(dft_section, "POISSON")
            CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
            ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
            print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
            CALL get_qs_kind_set(qs_kind_set, basis_rcut=ewald_rcut)
            CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat)
            CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
            CALL set_qs_env(qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw)
            CALL ewald_env_release(ewald_env)
            CALL ewald_pw_release(ewald_pw)
         END IF
      ELSEIF (dft_control%qs_control%method_id == do_method_xtb) THEN
         !   *** Read xTB parameter file ***
         xtb_control => dft_control%qs_control%xtb_control
         NULLIFY (ewald_env, ewald_pw)
         CALL get_qs_env(qs_env, nkind=nkind)
         DO ikind = 1, nkind
            qs_kind => qs_kind_set(ikind)
            ! Setup proper xTB parameters
            CPASSERT(.NOT. ASSOCIATED(qs_kind%xtb_parameter))
            CALL allocate_xtb_atom_param(qs_kind%xtb_parameter)
            ! Set default parameters
            CALL get_qs_kind(qs_kind, element_symbol=element_symbol)
            CALL xtb_parameters_init(qs_kind%xtb_parameter, element_symbol, &
                                     xtb_control%parameter_file_path, xtb_control%parameter_file_name, &
                                     para_env)
            ! Read specific parameters from input
            xtb_section => section_vals_get_subs_vals(qs_section, "xTB")
            CALL xtb_parameters_read(qs_kind%xtb_parameter, element_symbol, xtb_section)
            ! set dependent parameters
            CALL xtb_parameters_set(qs_kind%xtb_parameter)
            ! Generate basis set
            NULLIFY (tmp_basis_set)
            IF (qs_kind%xtb_parameter%z == 1) THEN
               ! special case hydrogen
               ngauss = xtb_control%h_sto_ng
            ELSE
               ngauss = xtb_control%sto_ng
            END IF
            CALL init_xtb_basis(qs_kind%xtb_parameter, tmp_basis_set, ngauss)
            CALL add_basis_set_to_container(qs_kind%basis_sets, tmp_basis_set, "ORB")
            ! potential
            zeff_correction = 0.0_dp
            CALL init_potential(qs_kind%all_potential, itype="BARE", &
                                zeff=qs_kind%xtb_parameter%zeff, zeff_correction=zeff_correction)
            CALL get_potential(qs_kind%all_potential, alpha_core_charge=alpha)
            ccore = qs_kind%xtb_parameter%zeff*SQRT((alpha/pi)**3)
            CALL set_potential(qs_kind%all_potential, ccore_charge=ccore)
            qs_kind%xtb_parameter%zeff = qs_kind%xtb_parameter%zeff - zeff_correction
         END DO
         !
         ! check for Ewald
         IF (xtb_control%do_ewald) THEN
            CALL ewald_env_create(ewald_env, para_env)
            poisson_section => section_vals_get_subs_vals(dft_section, "POISSON")
            CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
            ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
            print_section => section_vals_get_subs_vals(qs_env%input, "PRINT%GRID_INFORMATION")
            CALL read_ewald_section_tb(ewald_env, ewald_section, cell_ref%hmat)
            CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, print_section=print_section)
            CALL set_qs_env(qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw)
            CALL ewald_env_release(ewald_env)
            CALL ewald_pw_release(ewald_pw)
         END IF
      END IF

      ! lri or ri env initialization
      lri_section => section_vals_get_subs_vals(qs_section, "LRIGPW")
      IF (dft_control%qs_control%method_id == do_method_lrigpw .OR. &
          dft_control%qs_control%lri_optbas .OR. &
          dft_control%qs_control%method_id == do_method_rigpw) THEN
         CALL lri_env_init(lri_env, lri_section)
         CALL set_qs_env(qs_env, lri_env=lri_env)
      END IF

      !   *** Check basis and fill in missing parts ***
      CALL check_qs_kind_set(qs_kind_set, dft_control, subsys_section=subsys_section)

      !   *** Check that no all-electron potential is present if GPW or GAPW_XC
      CALL get_qs_kind_set(qs_kind_set, all_potential_present=all_potential_present)
      IF ((dft_control%qs_control%method_id == do_method_gpw) .OR. &
          (dft_control%qs_control%method_id == do_method_gapw_xc) .OR. &
          (dft_control%qs_control%method_id == do_method_ofgpw)) THEN
         IF (all_potential_present) THEN
            CPABORT("All-electron calculations with GPW, GAPW_XC, and OFGPW are not implemented")
         END IF
      END IF

      ! DFT+U
      CALL get_qs_kind_set(qs_kind_set, dft_plus_u_atom_present=dft_control%dft_plus_u)

      IF (dft_control%do_admm) THEN
         ! Check if ADMM basis is available, auto-generate if needed
         CALL get_qs_env(qs_env, nkind=nkind)
         DO ikind = 1, nkind
            NULLIFY (aux_fit_basis)
            qs_kind => qs_kind_set(ikind)
            CALL get_qs_kind(qs_kind, basis_set=aux_fit_basis, basis_type="AUX_FIT")
            IF (.NOT. (ASSOCIATED(aux_fit_basis))) THEN
               ! AUX_FIT basis set is not yet loaded
               CALL cp_warn(__LOCATION__, "Automatic Generation of AUX_FIT basis. "// &
                            "This is experimental code.")
               ! Generate a default basis
               CALL create_aux_fit_basis_set(aux_fit_basis, qs_kind, dft_control%auto_basis_aux_fit)
               CALL add_basis_set_to_container(qs_kind%basis_sets, aux_fit_basis, "AUX_FIT")
            END IF
         END DO
      END IF

      lribas = .FALSE.
      e1terms = .FALSE.
      IF (dft_control%qs_control%method_id == do_method_lrigpw) THEN
         lribas = .TRUE.
         CALL get_qs_env(qs_env, lri_env=lri_env)
         e1terms = lri_env%exact_1c_terms
      END IF
      IF (dft_control%qs_control%do_kg) THEN
         CALL section_vals_val_get(dft_section, "KG_METHOD%TNADD_METHOD", i_val=tnadd_method)
         IF (tnadd_method == kg_tnadd_embed_ri) lribas = .TRUE.
      END IF
      IF (lribas) THEN
         ! Check if LRI_AUX basis is available, auto-generate if needed
         CALL get_qs_env(qs_env, nkind=nkind)
         DO ikind = 1, nkind
            NULLIFY (lri_aux_basis)
            qs_kind => qs_kind_set(ikind)
            CALL get_qs_kind(qs_kind, basis_set=lri_aux_basis, basis_type="LRI_AUX")
            IF (.NOT. (ASSOCIATED(lri_aux_basis))) THEN
               ! LRI_AUX basis set is not yet loaded
               CALL cp_warn(__LOCATION__, "Automatic Generation of LRI_AUX basis. "// &
                            "This is experimental code.")
               ! Generate a default basis
               CALL create_lri_aux_basis_set(lri_aux_basis, qs_kind, dft_control%auto_basis_lri_aux, e1terms)
               CALL add_basis_set_to_container(qs_kind%basis_sets, lri_aux_basis, "LRI_AUX")
            END IF
         END DO
      END IF

      CALL section_vals_val_get(qs_env%input, "DFT%XC%HF%RI%_SECTION_PARAMETERS_", l_val=do_ri_hfx)
      IF (do_ri_hfx) THEN
         CALL get_qs_env(qs_env, nkind=nkind)
         DO ikind = 1, nkind
            NULLIFY (ri_hfx_basis)
            qs_kind => qs_kind_set(ikind)
            CALL get_qs_kind(qs_kind=qs_kind, basis_set=ri_hfx_basis, &
                             basis_type="RI_HFX")
            IF (.NOT. (ASSOCIATED(ri_hfx_basis))) THEN
               CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto)
               IF (dft_control%do_admm) THEN
                  CALL create_ri_aux_basis_set(ri_hfx_basis, qs_kind, dft_control%auto_basis_ri_hfx, basis_type="AUX_FIT")
               ELSE
                  CALL create_ri_aux_basis_set(ri_hfx_basis, qs_kind, dft_control%auto_basis_ri_hfx)
               ENDIF
               CALL add_basis_set_to_container(qs_kind%basis_sets, ri_hfx_basis, "RI_HFX")
            END IF
         END DO
      ENDIF

      IF (dft_control%qs_control%method_id == do_method_rigpw) THEN
         ! Check if RI_HXC basis is available, auto-generate if needed
         CALL get_qs_env(qs_env, nkind=nkind)
         DO ikind = 1, nkind
            NULLIFY (ri_hfx_basis)
            qs_kind => qs_kind_set(ikind)
            CALL get_qs_kind(qs_kind, basis_set=ri_hfx_basis, basis_type="RI_HXC")
            IF (.NOT. (ASSOCIATED(ri_hfx_basis))) THEN
               ! RI_HXC basis set is not yet loaded
               CALL cp_warn(__LOCATION__, "Automatic Generation of RI_HXC basis. "// &
                            "This is experimental code.")
               ! Generate a default basis
               CALL create_ri_aux_basis_set(ri_hfx_basis, qs_kind, dft_control%auto_basis_ri_hxc)
               CALL add_basis_set_to_container(qs_kind%basis_sets, ri_hfx_basis, "RI_HXC")
            END IF
         END DO
      END IF

      mp2_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION")
      CALL section_vals_get(mp2_section, explicit=mp2_present)
      IF (mp2_present) THEN

         ! basis should be sorted for imaginary time RPA/GW
         CALL section_vals_val_get(qs_env%input, "DFT%SORT_BASIS", i_val=sort_basis)
         CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%LOW_SCALING%_SECTION_PARAMETERS_", &
                                   l_val=do_wfc_im_time)

         IF (do_wfc_im_time .AND. sort_basis /= basis_sort_zet) THEN
            CALL cp_warn(__LOCATION__, &
                         "Low-scaling RPA requires SORT_BASIS EXP keyword (in DFT input section) for good performance")
         ENDIF

         ! Check if RI_AUX basis (for MP2/RPA) is given, auto-generate if not
         CALL mp2_env_create(qs_env%mp2_env)
         CALL get_qs_env(qs_env, mp2_env=mp2_env, nkind=nkind)
         CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_MP2%_SECTION_PARAMETERS_", l_val=do_ri_mp2)
         CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_SOS_MP2%_SECTION_PARAMETERS_", l_val=do_ri_sos_mp2)
         CALL section_vals_val_get(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%_SECTION_PARAMETERS_", l_val=do_ri_rpa)
         IF (do_ri_mp2 .OR. do_ri_sos_mp2 .OR. do_ri_rpa) THEN
            DO ikind = 1, nkind
               NULLIFY (ri_aux_basis_set)
               qs_kind => qs_kind_set(ikind)
               CALL get_qs_kind(qs_kind=qs_kind, basis_set=ri_aux_basis_set, &
                                basis_type="RI_AUX")
               IF (.NOT. (ASSOCIATED(ri_aux_basis_set))) THEN
                  ! RI_AUX basis set is not yet loaded
                  CALL cp_warn(__LOCATION__, "Automatic Generation of RI_AUX basis. "// &
                               "This is experimental code.")
                  ! Generate a default basis
                  CALL create_ri_aux_basis_set(ri_aux_basis_set, qs_kind, dft_control%auto_basis_ri_aux, basis_sort=sort_basis)
                  CALL add_basis_set_to_container(qs_kind%basis_sets, ri_aux_basis_set, "RI_AUX")
               END IF
            END DO
         END IF

      ENDIF

      IF (dft_control%do_xas_tdp_calculation) THEN
         ! Check if RI_XAS basis is given, auto-generate if not
         CALL get_qs_env(qs_env, nkind=nkind)
         DO ikind = 1, nkind
            NULLIFY (ri_xas_basis)
            qs_kind => qs_kind_set(ikind)
            CALL get_qs_kind(qs_kind, basis_Set=ri_xas_basis, basis_type="RI_XAS")
            IF (.NOT. ASSOCIATED(ri_xas_basis)) THEN
               CALL cp_warn(__LOCATION__, "Automatic Generation of RI_XAS basis. "// &
                            "This is experimental code.")
               ! Generate a default basis
               CALL create_ri_aux_basis_set(ri_xas_basis, qs_kind, dft_control%auto_basis_ri_xas)
               CALL add_basis_set_to_container(qs_kind%basis_sets, ri_xas_basis, "RI_XAS")
            END IF
         END DO
      END IF

      !   *** Initialize the spherical harmonics and ***
      !   *** the orbital transformation matrices    ***
      CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, maxlppl=maxlppl, maxlppnl=maxlppnl)

      lmax_sphere = dft_control%qs_control%gapw_control%lmax_sphere
      IF (lmax_sphere .LT. 0) THEN
         lmax_sphere = 2*maxlgto
         dft_control%qs_control%gapw_control%lmax_sphere = lmax_sphere
      END IF
      IF (dft_control%qs_control%method_id == do_method_lrigpw .OR. dft_control%qs_control%lri_optbas) THEN
         CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="LRI_AUX")
         !take maxlgto from lri basis if larger (usually)
         maxlgto = MAX(maxlgto, maxlgto_lri)
      ELSE IF (dft_control%qs_control%method_id == do_method_rigpw) THEN
         CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="RI_HXC")
         maxlgto = MAX(maxlgto, maxlgto_lri)
      END IF
      IF (dft_control%do_xas_tdp_calculation) THEN
         !done as a precaution
         CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_lri, basis_type="RI_XAS")
         maxlgto = MAX(maxlgto, maxlgto_lri)
      END IF
      maxl = MAX(2*maxlgto, maxlppl, maxlppnl, lmax_sphere) + 1

      CALL init_orbital_pointers(maxl)

      CALL init_spherical_harmonics(maxl, 0)

      !   *** Initialise the qs_kind_set ***
      CALL init_qs_kind_set(qs_kind_set)

      !   *** Initialise GAPW soft basis and projectors
      IF (dft_control%qs_control%method_id == do_method_gapw .OR. &
          dft_control%qs_control%method_id == do_method_gapw_xc) THEN
         qs_control => dft_control%qs_control
         CALL init_gapw_basis_set(qs_kind_set, qs_control, qs_env%input)
      ENDIF

      !   *** Initialize the pretabulation for the calculation of the   ***
      !   *** incomplete Gamma function F_n(t) after McMurchie-Davidson ***
      CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto)
      maxl = MAX(3*maxlgto + 1, 0)
      CALL init_md_ftable(maxl)

      !   *** Initialize the atomic interaction radii ***
      CALL init_interaction_radii(dft_control%qs_control, qs_kind_set)
      !
      IF (dft_control%qs_control%method_id == do_method_xtb) THEN
         ! cutoff radius
         DO ikind = 1, nkind
            qs_kind => qs_kind_set(ikind)
            CALL get_qs_kind(qs_kind, basis_set=tmp_basis_set)
            CALL get_gto_basis_set(tmp_basis_set, kind_radius=qs_kind%xtb_parameter%rcut)
         END DO
      END IF

      IF (.NOT. be_silent) THEN
         CALL write_pgf_orb_radii("orb", atomic_kind_set, qs_kind_set, subsys_section)
         CALL write_pgf_orb_radii("aux", atomic_kind_set, qs_kind_set, subsys_section)
         CALL write_pgf_orb_radii("lri", atomic_kind_set, qs_kind_set, subsys_section)
         CALL write_core_charge_radii(atomic_kind_set, qs_kind_set, subsys_section)
         CALL write_ppl_radii(atomic_kind_set, qs_kind_set, subsys_section)
         CALL write_ppnl_radii(atomic_kind_set, qs_kind_set, subsys_section)
         CALL write_paw_radii(atomic_kind_set, qs_kind_set, subsys_section)
      END IF

      !   *** Distribute molecules and atoms using the new data structures ***
      CALL distribute_molecules_1d(atomic_kind_set=atomic_kind_set, &
                                   particle_set=particle_set, &
                                   local_particles=local_particles, &
                                   molecule_kind_set=molecule_kind_set, &
                                   molecule_set=molecule_set, &
                                   local_molecules=local_molecules, &
                                   force_env_section=qs_env%input)

      !   *** SCF parameters ***
      CALL scf_c_create(scf_control)
      CALL scf_c_read_parameters(scf_control, dft_section)

      !   *** Allocate the data structure for Quickstep energies ***
      CALL allocate_qs_energy(energy)

      ! check for orthogonal basis
      has_unit_metric = .FALSE.
      IF (dft_control%qs_control%semi_empirical) THEN
         IF (dft_control%qs_control%se_control%orthogonal_basis) has_unit_metric = .TRUE.
      END IF
      IF (dft_control%qs_control%dftb) THEN
         IF (dft_control%qs_control%dftb_control%orthogonal_basis) has_unit_metric = .TRUE.
      END IF
      CALL set_qs_env(qs_env, has_unit_metric=has_unit_metric)

      !   *** Activate the interpolation ***
      CALL wfi_create(wf_history, &
                      interpolation_method_nr= &
                      dft_control%qs_control%wf_interpolation_method_nr, &
                      extrapolation_order=dft_control%qs_control%wf_extrapolation_order, &
                      has_unit_metric=has_unit_metric)

      !   *** Set the current Quickstep environment ***
      CALL set_qs_env(qs_env=qs_env, &
                      scf_control=scf_control, &
                      wf_history=wf_history)

      CALL qs_subsys_set(subsys, &
                         cell_ref=cell_ref, &
                         use_ref_cell=use_ref_cell, &
                         energy=energy, &
                         force=force)

      CALL get_qs_env(qs_env, ks_env=ks_env)
      CALL set_ks_env(ks_env, dft_control=dft_control)

      CALL qs_subsys_set(subsys, local_molecules=local_molecules, &
                         local_particles=local_particles, cell=cell)

      CALL distribution_1d_release(local_particles)
      CALL distribution_1d_release(local_molecules)
      CALL scf_c_release(scf_control)
      CALL wfi_release(wf_history)
      CALL dft_control_release(dft_control)

      CALL get_qs_env(qs_env=qs_env, &
                      atomic_kind_set=atomic_kind_set, &
                      dft_control=dft_control, &
                      scf_control=scf_control)

      ! decide what conditions need mo_derivs
      ! right now, this only appears to be OT
      IF (dft_control%qs_control%do_ls_scf .OR. &
          dft_control%qs_control%do_almo_scf) THEN
         CALL set_qs_env(qs_env=qs_env, requires_mo_derivs=.FALSE.)
      ELSE
         IF (scf_control%use_ot) THEN
            CALL set_qs_env(qs_env=qs_env, requires_mo_derivs=.TRUE.)
         ELSE
            CALL set_qs_env(qs_env=qs_env, requires_mo_derivs=.FALSE.)
         ENDIF
      ENDIF

      ! XXXXXXX this is backwards XXXXXXXX
      dft_control%smear = scf_control%smear%do_smear

      ! Periodic efield needs equal occupation and orbital gradients
      IF (.NOT. (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb)) THEN
         IF (dft_control%apply_period_efield) THEN
            CALL get_qs_env(qs_env=qs_env, requires_mo_derivs=orb_gradient)
            IF (.NOT. orb_gradient) THEN
               CALL cp_abort(__LOCATION__, "Periodic Efield needs orbital gradient and direct optimization."// &
                             " Use the OT optimization method.")
            END IF
            IF (dft_control%smear) THEN
               CALL cp_abort(__LOCATION__, "Periodic Efield needs equal occupation numbers."// &
                             " Smearing option is not possible.")
            END IF
         END IF
      END IF

      !   Initialize the GAPW local densities and potentials
      IF (dft_control%qs_control%method_id == do_method_gapw .OR. &
          dft_control%qs_control%method_id == do_method_gapw_xc) THEN
         !     *** Allocate and initialize the set of atomic densities ***
         NULLIFY (rho_atom_set)
         gapw_control => dft_control%qs_control%gapw_control
         CALL init_rho_atom(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
         CALL set_qs_env(qs_env=qs_env, rho_atom_set=rho_atom_set)
         IF (dft_control%qs_control%method_id /= do_method_gapw_xc) THEN
            CALL get_qs_env(qs_env=qs_env, local_rho_set=local_rho_set, natom=natom)
            !       *** Allocate and initialize the compensation density rho0 ***
            CALL init_rho0(local_rho_set, qs_env, gapw_control, .FALSE.)
            !       *** Allocate and Initialize the local coulomb term ***
            CALL init_coulomb_local(qs_env%hartree_local, natom)
         END IF
         ! NLCC
         CALL init_gapw_nlcc(qs_kind_set)
      ELSE IF (dft_control%qs_control%method_id == do_method_lrigpw) THEN
         ! allocate local ri environment
         ! nothing to do here?
      ELSE IF (dft_control%qs_control%method_id == do_method_rigpw) THEN
         ! allocate ri environment
         ! nothing to do here?
      ELSE IF (dft_control%qs_control%semi_empirical) THEN
         NULLIFY (se_store_int_env, se_nddo_mpole, se_nonbond_env)
         natom = SIZE(particle_set)
         se_section => section_vals_get_subs_vals(qs_section, "SE")
         se_control => dft_control%qs_control%se_control

         ! Make the cutoff radii choice a bit smarter
         CALL se_cutoff_compatible(se_control, se_section, cell, output_unit)

         SELECT CASE (dft_control%qs_control%method_id)
         CASE DEFAULT
         CASE (do_method_rm1, do_method_am1, do_method_mndo, do_method_pm3, &
               do_method_pm6, do_method_pm6fm, do_method_mndod, do_method_pnnl)
            ! Neighbor lists have to be MAX(interaction range, orbital range)
            ! set new kind radius
            CALL init_se_nlradius(se_control, atomic_kind_set, qs_kind_set, subsys_section)
         END SELECT
         ! Initialize to zero the max multipole to treat in the EWALD scheme..
         se_control%max_multipole = do_multipole_none
         ! check for Ewald
         IF (se_control%do_ewald .OR. se_control%do_ewald_gks) THEN
            CALL ewald_env_create(ewald_env, para_env)
            poisson_section => section_vals_get_subs_vals(dft_section, "POISSON")
            CALL ewald_env_set(ewald_env, poisson_section=poisson_section)
            ewald_section => section_vals_get_subs_vals(poisson_section, "EWALD")
            print_section => section_vals_get_subs_vals(qs_env%input, &
                                                        "PRINT%GRID_INFORMATION")
            CALL read_ewald_section(ewald_env, ewald_section)
            ! Create ewald grids
            CALL ewald_pw_create(ewald_pw, ewald_env, cell, cell_ref, &
                                 print_section=print_section)
            ! Initialize ewald grids
            CALL ewald_pw_grid_update(ewald_pw, ewald_env, cell%hmat)
            ! Setup the nonbond environment (real space part of Ewald)
            CALL ewald_env_get(ewald_env, rcut=ewald_rcut)
            ! Setup the maximum level of multipoles to be treated in the periodic SE scheme
            IF (se_control%do_ewald) THEN
               CALL ewald_env_get(ewald_env, max_multipole=se_control%max_multipole)
            ENDIF
            CALL section_vals_val_get(se_section, "NEIGHBOR_LISTS%VERLET_SKIN", &
                                      r_val=verlet_skin)
            CALL fist_nonbond_env_create(se_nonbond_env, atomic_kind_set, &
                                         do_nonbonded=.TRUE., verlet_skin=verlet_skin, ewald_rcut=ewald_rcut, &
                                         ei_scale14=0.0_dp, vdw_scale14=0.0_dp, shift_cutoff=.FALSE.)
            ! Create and Setup NDDO multipole environment
            CALL nddo_mpole_setup(se_nddo_mpole, natom)
            CALL set_qs_env(qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw, &
                            se_nonbond_env=se_nonbond_env, se_nddo_mpole=se_nddo_mpole)
            CALL ewald_env_release(ewald_env)
            CALL ewald_pw_release(ewald_pw)
            ! Handle the residual integral part 1/R^3
            CALL semi_empirical_expns3_setup(qs_kind_set, se_control, &
                                             dft_control%qs_control%method_id)
         END IF
         ! Taper function
         CALL se_taper_create(se_taper, se_control%integral_screening, se_control%do_ewald, &
                              se_control%taper_cou, se_control%range_cou, &
                              se_control%taper_exc, se_control%range_exc, &
                              se_control%taper_scr, se_control%range_scr, &
                              se_control%taper_lrc, se_control%range_lrc)
         CALL set_qs_env(qs_env, se_taper=se_taper)
         ! Store integral environment
         CALL semi_empirical_si_create(se_store_int_env, se_section)
         CALL set_qs_env(qs_env, se_store_int_env=se_store_int_env)
      ENDIF

      !   Initialize possible dispersion parameters
      IF (dft_control%qs_control%method_id == do_method_gpw .OR. &
          dft_control%qs_control%method_id == do_method_gapw .OR. &
          dft_control%qs_control%method_id == do_method_gapw_xc .OR. &
          dft_control%qs_control%method_id == do_method_lrigpw .OR. &
          dft_control%qs_control%method_id == do_method_rigpw .OR. &
          dft_control%qs_control%method_id == do_method_ofgpw) THEN
         ALLOCATE (dispersion_env)
         NULLIFY (xc_section)
         xc_section => section_vals_get_subs_vals(dft_section, "XC")
         CALL qs_dispersion_env_set(dispersion_env, xc_section)
         IF (dispersion_env%type == xc_vdw_fun_pairpot) THEN
            NULLIFY (pp_section)
            pp_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%PAIR_POTENTIAL")
            CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, pp_section, para_env)
         ELSE IF (dispersion_env%type == xc_vdw_fun_nonloc) THEN
            NULLIFY (nl_section)
            nl_section => section_vals_get_subs_vals(xc_section, "VDW_POTENTIAL%NON_LOCAL")
            CALL qs_dispersion_nonloc_init(dispersion_env, para_env)
         END IF
         CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
      ELSE IF (dft_control%qs_control%method_id == do_method_dftb) THEN
         ALLOCATE (dispersion_env)
         ! set general defaults
         dispersion_env%doabc = .FALSE.
         dispersion_env%c9cnst = .FALSE.
         dispersion_env%lrc = .FALSE.
         dispersion_env%srb = .FALSE.
         dispersion_env%verbose = .FALSE.
         NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, &
                  dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist, &
                  dispersion_env%d3_exclude_pair)
         NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, &
                  dispersion_env%d2y_dx2, dispersion_env%dftd_section)
         NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
         IF (dftb_control%dispersion .AND. dftb_control%dispersion_type == dispersion_d3) THEN
            dispersion_env%type = xc_vdw_fun_pairpot
            dispersion_env%pp_type = vdw_pairpot_dftd3
            dispersion_env%eps_cn = dftb_control%epscn
            dispersion_env%s6 = dftb_control%sd3(1)
            dispersion_env%sr6 = dftb_control%sd3(2)
            dispersion_env%s8 = dftb_control%sd3(3)
            dispersion_env%domol = .FALSE.
            dispersion_env%kgc8 = 0._dp
            dispersion_env%rc_disp = dftb_control%rcdisp
            dispersion_env%exp_pre = 0._dp
            dispersion_env%scaling = 0._dp
            dispersion_env%nd3_exclude_pair = 0
            dispersion_env%parameter_file_name = dftb_control%dispersion_parameter_file
            CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
         ELSE
            dispersion_env%type = xc_vdw_fun_none
         END IF
         CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
      ELSE IF (dft_control%qs_control%method_id == do_method_xtb) THEN
         ALLOCATE (dispersion_env)
         ! set general defaults
         dispersion_env%doabc = .FALSE.
         dispersion_env%c9cnst = .FALSE.
         dispersion_env%lrc = .FALSE.
         dispersion_env%srb = .FALSE.
         dispersion_env%verbose = .FALSE.
         NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, &
                  dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist, &
                  dispersion_env%d3_exclude_pair)
         NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, &
                  dispersion_env%d2y_dx2, dispersion_env%dftd_section)
         NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
         dispersion_env%type = xc_vdw_fun_pairpot
         dispersion_env%pp_type = vdw_pairpot_dftd3bj
         dispersion_env%eps_cn = xtb_control%epscn
         dispersion_env%s6 = xtb_control%s6
         dispersion_env%s8 = xtb_control%s8
         dispersion_env%a1 = xtb_control%a1
         dispersion_env%a2 = xtb_control%a2
         dispersion_env%domol = .FALSE.
         dispersion_env%kgc8 = 0._dp
         dispersion_env%rc_disp = xtb_control%rcdisp
         dispersion_env%exp_pre = 0._dp
         dispersion_env%scaling = 0._dp
         dispersion_env%nd3_exclude_pair = 0
         dispersion_env%parameter_file_name = xtb_control%dispersion_parameter_file
         CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
         CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
      ELSE IF (dft_control%qs_control%semi_empirical) THEN
         ALLOCATE (dispersion_env)
         ! set general defaults
         dispersion_env%doabc = .FALSE.
         dispersion_env%c9cnst = .FALSE.
         dispersion_env%lrc = .FALSE.
         dispersion_env%srb = .FALSE.
         dispersion_env%verbose = .FALSE.
         NULLIFY (dispersion_env%c6ab, dispersion_env%maxci, dispersion_env%r0ab, dispersion_env%rcov, &
                  dispersion_env%r2r4, dispersion_env%cn, dispersion_env%cnkind, dispersion_env%cnlist, &
                  dispersion_env%d3_exclude_pair)
         NULLIFY (dispersion_env%q_mesh, dispersion_env%kernel, dispersion_env%d2phi_dk2, &
                  dispersion_env%d2y_dx2, dispersion_env%dftd_section)
         NULLIFY (dispersion_env%sab_vdw, dispersion_env%sab_cn)
         IF (se_control%dispersion) THEN
            dispersion_env%type = xc_vdw_fun_pairpot
            dispersion_env%pp_type = vdw_pairpot_dftd3
            dispersion_env%eps_cn = se_control%epscn
            dispersion_env%s6 = se_control%sd3(1)
            dispersion_env%sr6 = se_control%sd3(2)
            dispersion_env%s8 = se_control%sd3(3)
            dispersion_env%domol = .FALSE.
            dispersion_env%kgc8 = 0._dp
            dispersion_env%rc_disp = se_control%rcdisp
            dispersion_env%exp_pre = 0._dp
            dispersion_env%scaling = 0._dp
            dispersion_env%nd3_exclude_pair = 0
            dispersion_env%parameter_file_name = se_control%dispersion_parameter_file
            CALL qs_dispersion_pairpot_init(atomic_kind_set, qs_kind_set, dispersion_env, para_env=para_env)
         ELSE
            dispersion_env%type = xc_vdw_fun_none
         END IF
         CALL set_qs_env(qs_env, dispersion_env=dispersion_env)
      END IF

      ! Initialize possible geomertical counterpoise correction potential
      IF (dft_control%qs_control%method_id == do_method_gpw .OR. &
          dft_control%qs_control%method_id == do_method_gapw .OR. &
          dft_control%qs_control%method_id == do_method_gapw_xc .OR. &
          dft_control%qs_control%method_id == do_method_lrigpw .OR. &
          dft_control%qs_control%method_id == do_method_rigpw .OR. &
          dft_control%qs_control%method_id == do_method_ofgpw) THEN
         ALLOCATE (gcp_env)
         NULLIFY (xc_section)
         xc_section => section_vals_get_subs_vals(dft_section, "XC")
         CALL qs_gcp_env_set(gcp_env, xc_section)
         CALL qs_gcp_init(qs_env, gcp_env)
         CALL set_qs_env(qs_env, gcp_env=gcp_env)
      END IF

      !   *** Allocate the MO data types ***
      CALL get_qs_kind_set(qs_kind_set, nsgf=n_ao, nelectron=nelectron)

      ! the total number of electrons
      nelectron = nelectron - dft_control%charge

      IF (dft_control%multiplicity == 0) THEN
         IF (MODULO(nelectron, 2) == 0) THEN
            dft_control%multiplicity = 1
         ELSE
            dft_control%multiplicity = 2
         END IF
      END IF

      multiplicity = dft_control%multiplicity

      IF ((dft_control%nspins < 1) .OR. (dft_control%nspins > 2)) THEN
         CPABORT("nspins should be 1 or 2 for the time being ...")
      END IF

      IF ((MODULO(nelectron, 2) /= 0) .AND. (dft_control%nspins == 1)) THEN
         IF (.NOT. dft_control%qs_control%ofgpw .AND. .NOT. dft_control%smear) THEN
            CPABORT("Use the LSD option for an odd number of electrons")
         END IF
      END IF

      ! The transition potential method to calculate XAS needs LSD
      IF (dft_control%do_xas_calculation) THEN
         IF (dft_control%nspins == 1) THEN
            CPABORT("Use the LSD option for XAS with transition potential")
         END IF
      END IF

      !   assigning the number of states per spin initial version, not yet very
      !   general. Should work for an even number of electrons and a single
      !   additional electron this set of options that requires full matrices,
      !   however, makes things a bit ugly right now.... we try to make a
      !   distinction between the number of electrons per spin and the number of
      !   MOs per spin this should allow the use of fractional occupations later
      !   on
      IF (dft_control%qs_control%ofgpw) THEN

         IF (dft_control%nspins == 1) THEN
            maxocc = nelectron
            nelectron_spin(1) = nelectron
            nelectron_spin(2) = 0
            n_mo(1) = 1
            n_mo(2) = 0
         ELSE
            IF (MODULO(nelectron + multiplicity - 1, 2) /= 0) THEN
               CPABORT("LSD: try to use a different multiplicity")
            END IF
            nelectron_spin(1) = (nelectron + multiplicity - 1)/2
            nelectron_spin(2) = (nelectron - multiplicity + 1)/2
            IF (nelectron_spin(1) < 0) THEN
               CPABORT("LSD: too few electrons for this multiplicity")
            END IF
            maxocc = MAXVAL(nelectron_spin)
            n_mo(1) = MIN(nelectron_spin(1), 1)
            n_mo(2) = MIN(nelectron_spin(2), 1)
         END IF

      ELSE

         IF (dft_control%nspins == 1) THEN
            maxocc = 2.0_dp
            nelectron_spin(1) = nelectron
            nelectron_spin(2) = 0
            IF (MODULO(nelectron, 2) == 0) THEN
               n_mo(1) = nelectron/2
            ELSE
               n_mo(1) = INT(nelectron/2._dp) + 1
            END IF
            n_mo(2) = 0
         ELSE
            maxocc = 1.0_dp

            ! The simplist spin distribution is written here. Special cases will
            ! need additional user input
            IF (MODULO(nelectron + multiplicity - 1, 2) /= 0) THEN
               CPABORT("LSD: try to use a different multiplicity")
            END IF

            nelectron_spin(1) = (nelectron + multiplicity - 1)/2
            nelectron_spin(2) = (nelectron - multiplicity + 1)/2

            IF (nelectron_spin(2) < 0) THEN
               CPABORT("LSD: too few electrons for this multiplicity")
            END IF

            n_mo(1) = nelectron_spin(1)
            n_mo(2) = nelectron_spin(2)

         END IF

      END IF

      ! Read the total_zeff_corr here [SGh]
      CALL get_qs_kind_set(qs_kind_set, total_zeff_corr=total_zeff_corr)
      ! store it in qs_env
      qs_env%total_zeff_corr = total_zeff_corr

      ! store the number of electrons once an for all
      CALL qs_subsys_set(subsys, &
                         nelectron_total=nelectron, &
                         nelectron_spin=nelectron_spin)

      ! Check and set number of added (unoccupied) MOs
      IF (dft_control%nspins == 2) THEN
         IF (scf_control%added_mos(2) > 0) THEN
            n_mo_add = scf_control%added_mos(2)
         ELSE
            n_mo_add = scf_control%added_mos(1)
         END IF
         IF (n_mo_add > n_ao - n_mo(2)) &
            CPWARN("More added MOs requested for beta spin than available.")
         scf_control%added_mos(2) = MIN(n_mo_add, n_ao - n_mo(2))
         n_mo(2) = n_mo(2) + scf_control%added_mos(2)
      END IF

      ! proceed alpha orbitals after the beta orbitals; this is essential to avoid
      ! reduction in the number of available unoccupied molecular orbitals.
      ! E.g. n_ao = 10, nelectrons = 10, multiplicity = 3 implies n_mo(1) = 6, n_mo(2) = 4;
      ! added_mos(1:2) = (6,undef) should increase the number of molecular orbitals as
      ! n_mo(1) = min(n_ao, n_mo(1) + added_mos(1)) = 10, n_mo(2) = 10.
      ! However, if we try to proceed alpha orbitals first, this leads us n_mo(1:2) = (10,8)
      ! due to the following assignment instruction above:
      !   IF (scf_control%added_mos(2) > 0) THEN ... ELSE; n_mo_add = scf_control%added_mos(1); END IF
      IF (scf_control%added_mos(1) > n_ao - n_mo(1)) &
         CALL cp_warn(__LOCATION__, &
                      "More added MOs requested than available. "// &
                      "The full set of unoccupied MOs will be used.")
      scf_control%added_mos(1) = MIN(scf_control%added_mos(1), n_ao - n_mo(1))
      n_mo(1) = n_mo(1) + scf_control%added_mos(1)

      IF (dft_control%nspins == 2) THEN
         IF (n_mo(2) > n_mo(1)) &
            CALL cp_warn(__LOCATION__, &
                         "More beta than alpha MOs requested. "// &
                         "The number of beta MOs will be reduced to the number alpha MOs.")
         n_mo(2) = MIN(n_mo(1), n_mo(2))
         CPASSERT(n_mo(1) >= nelectron_spin(1))
         CPASSERT(n_mo(2) >= nelectron_spin(2))
      END IF

      ! kpoints
      CALL get_qs_env(qs_env=qs_env, do_kpoints=do_kpoints)
      IF (do_kpoints .AND. dft_control%nspins == 2) THEN
         ! we need equal number of calculated states
         IF (n_mo(2) /= n_mo(1)) &
            CALL cp_warn(__LOCATION__, &
                         "Kpoints: Different number of MOs requested. "// &
                         "The number of beta MOs will be set to the number alpha MOs.")
         n_mo(2) = n_mo(1)
         CPASSERT(n_mo(1) >= nelectron_spin(1))
         CPASSERT(n_mo(2) >= nelectron_spin(2))
      END IF

      ! Compatibility checks for smearing
      IF (scf_control%smear%do_smear) THEN
         IF (scf_control%added_mos(1) == 0) THEN
            CPABORT("Extra MOs (ADDED_MOS) are required for smearing")
         END IF
      END IF

      !   *** Some options require that all MOs are computed ... ***
      IF (BTEST(cp_print_key_should_output(logger%iter_info, dft_section, &
                                           "PRINT%MO/CARTESIAN"), &
                cp_p_file) .OR. &
          (scf_control%level_shift /= 0.0_dp) .OR. &
          (scf_control%diagonalization%eps_jacobi /= 0.0_dp) .OR. &
          (dft_control%roks .AND. (.NOT. scf_control%use_ot))) THEN
         n_mo(:) = n_ao
      END IF

      ! Compatibility checks for ROKS
      IF (dft_control%roks .AND. (.NOT. scf_control%use_ot)) THEN
         IF (scf_control%roks_scheme == general_roks) &
            CPWARN("General ROKS scheme is not yet tested!")

         IF (scf_control%smear%do_smear) THEN
            CALL cp_abort(__LOCATION__, &
                          "The options ROKS and SMEAR are not compatible. "// &
                          "Try UKS instead of ROKS")
         END IF
      END IF

      ! in principle the restricted calculation could be performed
      ! using just one set of MOs and special casing most of the code
      ! right now we'll just take care of what is effectively an additional constraint
      ! at as few places as possible, just duplicating the beta orbitals
      IF (dft_control%restricted .AND. (output_unit > 0)) THEN
         ! it is really not yet tested till the end ! Joost
         WRITE (output_unit, *) ""
         WRITE (output_unit, *) " **************************************"
         WRITE (output_unit, *) " restricted calculation cutting corners"
         WRITE (output_unit, *) " experimental feature, check code      "
         WRITE (output_unit, *) " **************************************"
      ENDIF

      ! no point in allocating these things here ?
      IF (dft_control%qs_control%do_ls_scf) THEN
         NULLIFY (mos)
      ELSE
         ALLOCATE (mos(dft_control%nspins))
         DO ispin = 1, dft_control%nspins
            NULLIFY (mos(ispin)%mo_set)
            CALL allocate_mo_set(mo_set=mos(ispin)%mo_set, &
                                 nao=n_ao, &
                                 nmo=n_mo(ispin), &
                                 nelectron=nelectron_spin(ispin), &
                                 n_el_f=REAL(nelectron_spin(ispin), dp), &
                                 maxocc=maxocc, &
                                 flexible_electron_count=dft_control%relax_multiplicity)
         END DO
      END IF

      CALL set_qs_env(qs_env, mos=mos)

      ! If we use auxiliary density matrix methods , set mo_set_aux_fit
      IF (dft_control%do_admm) THEN
         !
         CALL get_qs_kind_set(qs_kind_set, nelectron=nelectron, nsgf=n_ao_aux_fit, &
                              basis_type="AUX_FIT")
         ALLOCATE (mos_aux_fit(dft_control%nspins))

         DO ispin = 1, dft_control%nspins
            NULLIFY (mos_aux_fit(ispin)%mo_set)
            CALL allocate_mo_set(mo_set=mos_aux_fit(ispin)%mo_set, &
                                 nao=n_ao_aux_fit, &
                                 nmo=n_mo(ispin), &
                                 nelectron=nelectron_spin(ispin), &
                                 n_el_f=REAL(nelectron_spin(ispin), dp), &
                                 maxocc=maxocc, &
                                 flexible_electron_count=dft_control%relax_multiplicity)
         END DO
         CALL set_qs_env(qs_env, mos_aux_fit=mos_aux_fit)
      END IF

! allocate mos when switch_surf_dip is triggered [SGh]
      IF (dft_control%switch_surf_dip) THEN
         ALLOCATE (mos_last_converged(dft_control%nspins))
         DO ispin = 1, dft_control%nspins
            NULLIFY (mos_last_converged(ispin)%mo_set)
            CALL allocate_mo_set(mo_set=mos_last_converged(ispin)%mo_set, &
                                 nao=n_ao, &
                                 nmo=n_mo(ispin), &
                                 nelectron=nelectron_spin(ispin), &
                                 n_el_f=REAL(nelectron_spin(ispin), dp), &
                                 maxocc=maxocc, &
                                 flexible_electron_count=dft_control%relax_multiplicity)
         END DO
         CALL set_qs_env(qs_env, mos_last_converged=mos_last_converged)
      ENDIF

      IF (.NOT. be_silent) THEN
         ! Print the DFT control parameters
         CALL write_dft_control(dft_control, dft_section)

         ! Print the vdW control parameters
         IF (dft_control%qs_control%method_id == do_method_gpw .OR. &
             dft_control%qs_control%method_id == do_method_gapw .OR. &
             dft_control%qs_control%method_id == do_method_gapw_xc .OR. &
             dft_control%qs_control%method_id == do_method_lrigpw .OR. &
             dft_control%qs_control%method_id == do_method_rigpw .OR. &
             dft_control%qs_control%method_id == do_method_dftb .OR. &
             dft_control%qs_control%method_id == do_method_xtb .OR. &
             dft_control%qs_control%method_id == do_method_ofgpw) THEN
            CALL get_qs_env(qs_env, dispersion_env=dispersion_env)
            CALL qs_write_dispersion(qs_env, dispersion_env)
         END IF

         ! Print the Quickstep control parameters
         CALL write_qs_control(dft_control%qs_control, dft_section)

         ! Print XES/XAS control parameters
         IF (dft_control%do_xas_calculation) THEN
            CALL cite_reference(Iannuzzi2007)
            !CALL write_xas_control(dft_control%xas_control,dft_section)
         END IF

         ! Print the unnormalized basis set information (input data)
         CALL write_gto_basis_sets(qs_kind_set, subsys_section)

         ! Print the atomic kind set
         CALL write_qs_kind_set(qs_kind_set, subsys_section)

         ! Print the molecule kind set
         CALL write_molecule_kind_set(molecule_kind_set, subsys_section)

         ! Print the total number of kinds, atoms, basis functions etc.
         CALL write_total_numbers(qs_kind_set, particle_set, qs_env%input)

         ! Print the atomic coordinates
         CALL write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label="QUICKSTEP")

         ! Print the interatomic distances
         CALL write_particle_distances(particle_set, cell, subsys_section)

         ! Print the requested structure data
         CALL write_structure_data(particle_set, cell, subsys_section)

         ! Print symmetry information
         CALL write_symmetry(particle_set, cell, subsys_section)

         ! Print the SCF parameters
         IF ((.NOT. dft_control%qs_control%do_ls_scf) .AND. &
             (.NOT. dft_control%qs_control%do_almo_scf)) THEN
            CALL scf_c_write_parameters(scf_control, dft_section)
         ENDIF
      ENDIF

      ! Sets up pw_env, qs_charges, mpools ...
      CALL qs_env_setup(qs_env)

      ! Allocate and initialise rho0 soft on the global grid
      IF (dft_control%qs_control%method_id == do_method_gapw) THEN
         CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho0_mpole=rho0_mpole)
         CALL rho0_s_grid_create(pw_env, rho0_mpole)
      END IF

      IF (output_unit > 0) CALL m_flush(output_unit)
      CALL timestop(handle)

   END SUBROUTINE qs_init_subsys

! **************************************************************************************************
!> \brief Write the total number of kinds, atoms, etc. to the logical unit
!>      number lunit.
!> \param qs_kind_set ...
!> \param particle_set ...
!> \param force_env_section ...
!> \author Creation (06.10.2000)
! **************************************************************************************************
   SUBROUTINE write_total_numbers(qs_kind_set, particle_set, force_env_section)

      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      TYPE(section_vals_type), POINTER                   :: force_env_section

      INTEGER                                            :: maxlgto, maxlppl, maxlppnl, natom, ncgf, &
                                                            nkind, npgf, nset, nsgf, nshell, &
                                                            output_unit
      TYPE(cp_logger_type), POINTER                      :: logger

      NULLIFY (logger)
      logger => cp_get_default_logger()
      output_unit = cp_print_key_unit_nr(logger, force_env_section, "PRINT%TOTAL_NUMBERS", &
                                         extension=".Log")

      IF (output_unit > 0) THEN
         natom = SIZE(particle_set)
         nkind = SIZE(qs_kind_set)

         CALL get_qs_kind_set(qs_kind_set, &
                              maxlgto=maxlgto, &
                              ncgf=ncgf, &
                              npgf=npgf, &
                              nset=nset, &
                              nsgf=nsgf, &
                              nshell=nshell, &
                              maxlppl=maxlppl, &
                              maxlppnl=maxlppnl)

         WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") &
            "TOTAL NUMBERS AND MAXIMUM NUMBERS"

         IF (nset + npgf + ncgf > 0) THEN
            WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T71,I10))") &
               "Total number of", &
               "- Atomic kinds:                  ", nkind, &
               "- Atoms:                         ", natom, &
               "- Shell sets:                    ", nset, &
               "- Shells:                        ", nshell, &
               "- Primitive Cartesian functions: ", npgf, &
               "- Cartesian basis functions:     ", ncgf, &
               "- Spherical basis functions:     ", nsgf
         ELSE IF (nshell + nsgf > 0) THEN
            WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T71,I10))") &
               "Total number of", &
               "- Atomic kinds:                  ", nkind, &
               "- Atoms:                         ", natom, &
               "- Shells:                        ", nshell, &
               "- Spherical basis functions:     ", nsgf
         ELSE
            WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T71,I10))") &
               "Total number of", &
               "- Atomic kinds:                  ", nkind, &
               "- Atoms:                         ", natom
         END IF

         IF ((maxlppl > -1) .AND. (maxlppnl > -1)) THEN
            WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T75,I6))") &
               "Maximum angular momentum of the", &
               "- Orbital basis functions:                   ", maxlgto, &
               "- Local part of the GTH pseudopotential:     ", maxlppl, &
               "- Non-local part of the GTH pseudopotential: ", maxlppnl
         ELSEIF (maxlppl > -1) THEN
            WRITE (UNIT=output_unit, FMT="(/,T3,A,(T30,A,T75,I6))") &
               "Maximum angular momentum of the", &
               "- Orbital basis functions:                   ", maxlgto, &
               "- Local part of the GTH pseudopotential:     ", maxlppl
         ELSE
            WRITE (UNIT=output_unit, FMT="(/,T3,A,T75,I6)") &
               "Maximum angular momentum of the orbital basis functions: ", maxlgto
         END IF

         ! LRI_AUX BASIS
         CALL get_qs_kind_set(qs_kind_set, &
                              maxlgto=maxlgto, &
                              ncgf=ncgf, &
                              npgf=npgf, &
                              nset=nset, &
                              nsgf=nsgf, &
                              nshell=nshell, &
                              basis_type="LRI_AUX")
         IF (nset + npgf + ncgf > 0) THEN
            WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,(T30,A,T71,I10))") &
               "LRI_AUX Basis: ", &
               "Total number of", &
               "- Shell sets:                    ", nset, &
               "- Shells:                        ", nshell, &
               "- Primitive Cartesian functions: ", npgf, &
               "- Cartesian basis functions:     ", ncgf, &
               "- Spherical basis functions:     ", nsgf
            WRITE (UNIT=output_unit, FMT="(T30,A,T75,I6)") &
               "  Maximum angular momentum ", maxlgto
         END IF

         ! RI_HXC BASIS
         CALL get_qs_kind_set(qs_kind_set, &
                              maxlgto=maxlgto, &
                              ncgf=ncgf, &
                              npgf=npgf, &
                              nset=nset, &
                              nsgf=nsgf, &
                              nshell=nshell, &
                              basis_type="RI_HXC")
         IF (nset + npgf + ncgf > 0) THEN
            WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,(T30,A,T71,I10))") &
               "RI_HXC Basis: ", &
               "Total number of", &
               "- Shell sets:                    ", nset, &
               "- Shells:                        ", nshell, &
               "- Primitive Cartesian functions: ", npgf, &
               "- Cartesian basis functions:     ", ncgf, &
               "- Spherical basis functions:     ", nsgf
            WRITE (UNIT=output_unit, FMT="(T30,A,T75,I6)") &
               "  Maximum angular momentum ", maxlgto
         END IF

         ! AUX_FIT BASIS
         CALL get_qs_kind_set(qs_kind_set, &
                              maxlgto=maxlgto, &
                              ncgf=ncgf, &
                              npgf=npgf, &
                              nset=nset, &
                              nsgf=nsgf, &
                              nshell=nshell, &
                              basis_type="AUX_FIT")
         IF (nset + npgf + ncgf > 0) THEN
            WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,(T30,A,T71,I10))") &
               "AUX_FIT ADMM-Basis: ", &
               "Total number of", &
               "- Shell sets:                    ", nset, &
               "- Shells:                        ", nshell, &
               "- Primitive Cartesian functions: ", npgf, &
               "- Cartesian basis functions:     ", ncgf, &
               "- Spherical basis functions:     ", nsgf
            WRITE (UNIT=output_unit, FMT="(T30,A,T75,I6)") &
               "  Maximum angular momentum ", maxlgto
         END IF

      END IF
      CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
                                        "PRINT%TOTAL_NUMBERS")

   END SUBROUTINE write_total_numbers

END MODULE qs_environment
