Initialization

  recursive function multisystem_basic_constructor(namespace, factory) result(system)
    type(namespace_t),            intent(in) :: namespace
    class(system_factory_abst_t), intent(in) :: factory
    class(multisystem_basic_t),   pointer    :: system

    PUSH_SUB(multisystem_basic_constructor)

    SAFE_ALLOCATE(system)

    call multisystem_basic_init(system, namespace, factory)

    POP_SUB(multisystem_basic_constructor)
  end function multisystem_basic_constructor

  recursive subroutine multisystem_basic_init(this, namespace, factory)
    class(multisystem_t),      intent(inout) :: this
    type(namespace_t),            intent(in) :: namespace
    class(system_factory_abst_t), intent(in) :: factory

    PUSH_SUB(multisystem_basic_init)

    call multisystem_init(this, namespace, factory)

    POP_SUB(multisystem_basic_init)
  end subroutine multisystem_basic_init

  recursive subroutine multisystem_init(this, namespace, factory)
    class(multisystem_t),      intent(inout) :: this
    type(namespace_t),            intent(in) :: namespace
    class(system_factory_abst_t), intent(in) :: factory

    integer :: isys, system_type, ic
    character(len=128) :: system_name
    type(block_t) :: blk

    PUSH_SUB(multisystem_init)

    this%namespace = namespace

    if (parse_block(this%namespace, factory%block_name(), blk) == 0) then

      do isys = 1, parse_block_n(blk)
        ! Parse system name and type
        call parse_block_string(blk, isys - 1, 0, system_name)
        if (len_trim(system_name) == 0) then
          call messages_input_error(this%namespace, factory%block_name(), 'All systems must have a name')
        end if
        do ic = 1, len(parser_varname_excluded_characters)
          if (index(trim(system_name), parser_varname_excluded_characters(ic:ic)) /= 0) then
            call messages_input_error(this%namespace, factory%block_name(), &
              'Illegal character "' // parser_varname_excluded_characters(ic:ic) // '" in system name', row=isys-1, column=0)
          end if
        end do
        call parse_block_integer(blk, isys - 1, 1, system_type)

        call multisystem_create_system(this, system_name, system_type, isys, factory)
      end do
      call parse_block_end(blk)
    else
      message(1) = "Input error while reading block "//trim(this%namespace%get())//"."//trim(factory%block_name())
      call messages_fatal(1, namespace=this%namespace)
    end if

    POP_SUB(multisystem_init)
  end subroutine multisystem_init

  recursive subroutine multisystem_create_system(this, system_name, system_type, isys, factory)
    class(multisystem_t),      intent(inout) :: this
    character(len=128),           intent(in) :: system_name
    integer,                      intent(in) :: system_type
    integer,                      intent(in) :: isys
    class(system_factory_abst_t), intent(in) :: factory

    type(system_iterator_t) :: iter
    class(system_t), pointer :: sys, other

    PUSH_SUB(multisystem_create_system)

    ! Create folder to store system files.
    ! Needs to be done before creating the system as this in turn might create subfolders.
    call io_mkdir(system_name, namespace=this%namespace)

    ! Create system
    sys => factory%create(this%namespace, system_name, system_type)
    if (.not. associated(sys)) then
      call messages_input_error(this%namespace, factory%block_name(), 'Unknown system type.')
    end if

    ! Check that the system is unique
    call iter%start(this%list)
    do while (iter%has_next())
      other => iter%get_next()
      if (sys%namespace == other%namespace) then
        call messages_input_error(this%namespace, factory%block_name(), 'Duplicated system in multisystem', &
          row=isys-1, column=0)
      end if
    end do

    ! Add system to list of systems
    call this%list%add(sys)

    POP_SUB(multisystem_create_system)
  end subroutine multisystem_create_system

  recursive subroutine multisystem_init_propagator(this)
    class(multisystem_t),      intent(inout) :: this

    type(system_iterator_t) :: iter
    class(system_t), pointer :: system
    type(interaction_iterator_t) :: inter_iter
    class(interaction_t), pointer :: interaction

    PUSH_SUB(multisystem_init_propagator)

    ! Now initialized the propagators of the subsystems
    call iter%start(this%list)
    do while (iter%has_next())
      system => iter%get_next()
      call system%init_propagator()
    end do

    ! Initialize the propagator of the multisystem. By default the
    ! multisystem itself and its own quantities are kept unchaged
    ! by using the static propagator. However, the subsystems are allowed to have
    ! their own propagators and those do not have to be static.
    ! Needs to be done after initializing the subsystems propagators,
    ! as we use the smallest dt of the subsystems.
    this%prop => propagator_static_t(this%smallest_algo_dt())
    this%interaction_timing = OPTION__INTERACTIONTIMING__TIMING_EXACT
    call this%prop%rewind()

    ! Initialize propagator clock
    this%prop%clock = clock_t(time_step=this%prop%dt/this%prop%algo_steps)

    ! Initialize system clock
    this%clock = clock_t(time_step=this%prop%dt)

    ! Interaction clocks
    call inter_iter%start(this%interactions)
    do while (inter_iter%has_next())
      interaction => inter_iter%get_next()
      interaction%clock = this%prop%clock - CLOCK_TICK
    end do

    ! Required quantities clocks
    where (this%quantities%required)
      this%quantities%clock = this%prop%clock
    end where

    POP_SUB(multisystem_init_propagator)
  end subroutine multisystem_init_propagator

  subroutine system_init_propagator(this)
    class(system_t),      intent(inout) :: this

    integer :: prop
    FLOAT :: dt
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction

    PUSH_SUB(system_init_propagator)

    call messages_experimental('Multisystem propagator framework')

    !%Variable TDSystemPropagator
    !%Type integer
    !%Default static
    !%Section Time-Dependent::Propagation
    !%Description
    !% A variable to set the propagator in the multisystem framework.
    !% This is a temporary solution, and should be replaced by the
    !% TDPropagator variable.
    !%Option static 0
    !% (Experimental) Do not propagate the system in time.
    !%Option verlet 1
    !% (Experimental) Verlet propagator.
    !%Option beeman 2
    !% (Experimental) Beeman propagator without predictor-corrector.
    !%Option beeman_scf 3
    !% (Experimental) Beeman propagator with predictor-corrector scheme.
    !%Option exp_mid 4
    !% (Experimental) Exponential midpoint propagator without predictor-corrector.
    !%Option exp_mid_scf 5
    !% (Experimental) Exponential midpoint propagator with predictor-corrector scheme.
    !%End
    call parse_variable(this%namespace, 'TDSystemPropagator', PROP_STATIC, prop)
    if (.not.varinfo_valid_option('TDSystemPropagator', prop)) call messages_input_error(this%namespace, 'TDSystemPropagator')
    call messages_print_var_option(stdout, 'TDSystemPropagator', prop)

    ! This variable is also defined (and properly documented) in td/td.F90.
    ! This is temporary, until all the propagators are moved to the new framework.
    call parse_variable(this%namespace, 'TDTimeStep', CNST(10.0), dt)
    if (dt <= M_ZERO) then
      call messages_input_error(this%namespace, 'TDTimeStep', "must be greater than zero")
    end if
    call messages_print_var_option(stdout, 'TDSystemPropagator', prop)

    select case (prop)
    case (PROP_STATIC)
      this%prop => propagator_static_t(dt)
    case (PROP_VERLET)
      this%prop => propagator_verlet_t(dt)
    case (PROP_BEEMAN)
      this%prop => propagator_beeman_t(dt, predictor_corrector=.false.)
    case (PROP_BEEMAN_SCF)
      this%prop => propagator_beeman_t(dt, predictor_corrector=.true.)
    case (PROP_EXPMID)
      this%prop => propagator_exp_mid_t(dt, predictor_corrector=.false.)
    case (PROP_EXPMID_SCF)
      this%prop => propagator_exp_mid_t(dt, predictor_corrector=.true.)
    case default
      call messages_input_error(this%namespace, 'TDSystemPropagator')
    end select

    call this%prop%rewind()

    ! Initialize propagator clock
    this%prop%clock = clock_t(time_step=this%prop%dt/this%prop%algo_steps)

    ! Initialize system clock
    this%clock = clock_t(time_step=this%prop%dt)

    ! Interaction clocks
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      interaction%clock = this%prop%clock - CLOCK_TICK
    end do

    ! Required quantities clocks
    where (this%quantities%required)
      this%quantities%clock = this%prop%clock
    end where

    !%Variable InteractionTiming
    !%Type integer
    !%Default timing_exact
    !%Section Time-Dependent::Propagation
    !%Description
    !% A parameter to determine if interactions should use the quantities
    !% at the exact time or if retardation is allowed.
    !%Option timing_exact 1
    !% Only allow interactions at exactly the same times
    !%Option timing_retarded 2
    !% Allow retarded interactions
    !%End
    call parse_variable(this%namespace, 'InteractionTiming', &
      OPTION__INTERACTIONTIMING__TIMING_EXACT, &
      this%interaction_timing)
    if (.not.varinfo_valid_option('InteractionTiming', this%interaction_timing)) then
      call messages_input_error(this%namespace, 'InteractionTiming')
    end if
    call messages_print_var_option(stdout, 'InteractionTiming', &
      this%interaction_timing)

    POP_SUB(system_init_propagator)
  end subroutine system_init_propagator