system.F90

Go to the documentation of this file.

!! Copyright (C) 2019 N. Tancogne-Dejean
!! Copyright (C) 2020 M. Oliveira, Heiko Appel
!! Copyright (C) 2021 S. Ohlmann
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module system_oct_m
  use algorithm_oct_m
  use algorithm_factory_oct_m
  use clock_oct_m
  use debug_oct_m
  use ghost_interaction_oct_m
  use global_oct_m
  use interaction_oct_m
  use interaction_partner_oct_m
  use interaction_with_partner_oct_m
  use messages_oct_m
  use mpi_oct_m
  use namespace_oct_m
  use multisystem_debug_oct_m
  use linked_list_oct_m
  use parser_oct_m
  use profiling_oct_m
  use quantity_oct_m
  use space_oct_m
  use unit_oct_m
  use unit_system_oct_m
  use varinfo_oct_m
  implicit none
 
  private
  public ::                                    &
    system_t,                                  &
    system_execute_algorithm,                  &
    system_init_parallelization,               &
    system_init_algorithm,                     &
    system_init_clocks,                        &
    system_reset_clocks,                       &
    system_propagation_start,                  &
    system_propagation_finish,                 &
    system_restart_read,                       &
    system_restart_write,                      &
    system_update_potential_energy,            &
    system_update_total_energy,                &
    system_end,                                &
    system_list_t,                             &
    system_iterator_t
 
  type :: barrier_t
    logical :: active
    float :: target_time
  end type barrier_t
 
  integer, parameter, public :: &
    NUMBER_BARRIERS = 1,        &
    barrier_restart = 1
 
  type, extends(interaction_partner_t), abstract :: system_t
    private
    class(algorithm_t),  pointer, public :: algo => null()
 
    integer, public :: accumulated_loop_ticks
 
    integer, public :: interaction_timing
 
    type(integer_list_t), public :: supported_interactions
    type(interaction_list_t), public :: interactions
 
    type(mpi_grp_t), public :: grp
 
    type(barrier_t) :: barrier(NUMBER_BARRIERS)
    float, public :: kinetic_energy     
    float, public :: potential_energy   
    float, public :: internal_energy    
    float, public :: total_energy       
 
  contains
    procedure :: execute_algorithm =>  system_execute_algorithm 
    procedure :: reset_clocks => system_reset_clocks 
    procedure :: update_exposed_quantities => system_update_exposed_quantities 
    procedure :: init_algorithm => system_init_algorithm 
    procedure :: algorithm_finished => system_algorithm_finished 
    procedure :: init_clocks => system_init_clocks 
    procedure :: init_all_interactions => system_init_all_interactions 
    procedure :: init_parallelization => system_init_parallelization 
    procedure :: update_interactions => system_update_interactions 
    procedure :: update_interactions_start => system_update_interactions_start 
    procedure :: update_interactions_finish => system_update_interactions_finish 
    procedure :: propagation_start => system_propagation_start 
    procedure :: propagation_finish => system_propagation_finish 
    procedure :: iteration_info => system_iteration_info 
    procedure :: restart_write => system_restart_write 
    procedure :: restart_read => system_restart_read 
    procedure :: output_start => system_output_start 
    procedure :: output_write => system_output_write 
    procedure :: output_finish => system_output_finish 
    procedure :: process_is_slave => system_process_is_slave 
    procedure :: start_barrier => system_start_barrier 
    procedure :: end_barrier => system_end_barrier 
    procedure :: arrived_at_barrier => system_arrived_at_barrier 
    procedure :: arrived_at_any_barrier => system_arrived_at_any_barrier 
    procedure :: update_potential_energy => system_update_potential_energy 
    procedure :: update_internal_energy => system_update_internal_energy 
    procedure :: update_total_energy => system_update_total_energy 
    procedure(system_init_interaction),          deferred :: init_interaction
    procedure(system_initial_conditions),        deferred :: initial_conditions
    procedure(system_do_algorithmic_operation),  deferred :: do_algorithmic_operation
    procedure(system_is_tolerance_reached),      deferred :: is_tolerance_reached
    procedure(system_update_quantity),           deferred :: update_quantity
    procedure(system_update_exposed_quantity),   deferred :: update_exposed_quantity
    procedure(system_restart_write_data),        deferred :: restart_write_data
    procedure(system_restart_read_data),         deferred :: restart_read_data
    procedure(system_update_kinetic_energy),     deferred :: update_kinetic_energy
  end type system_t
 
  abstract interface
 
    ! ---------------------------------------------------------
    subroutine system_init_interaction(this, interaction)
      import system_t
      import interaction_t
      class(system_t), target, intent(inout) :: this
      class(interaction_t),    intent(inout) :: interaction
    end subroutine system_init_interaction
 
    ! ---------------------------------------------------------
    subroutine system_initial_conditions(this)
      import system_t
      class(system_t), intent(inout) :: this
    end subroutine system_initial_conditions
 
    ! ---------------------------------------------------------
    logical function system_do_algorithmic_operation(this, operation, updated_quantities) result(done)
      import system_t
      import algorithmic_operation_t
      class(system_t),                intent(inout) :: this
      class(algorithmic_operation_t), intent(in)    :: operation
      integer, allocatable,           intent(out)   :: updated_quantities(:)
    end function system_do_algorithmic_operation
 
    ! ---------------------------------------------------------
    logical function system_is_tolerance_reached(this, tol)
      import system_t
      class(system_t), intent(in) :: this
      float,           intent(in) :: tol
    end function system_is_tolerance_reached
 
    ! ---------------------------------------------------------
    subroutine system_store_current_status(this)
      import system_t
      class(system_t), intent(inout) :: this
    end subroutine system_store_current_status
 
    ! ---------------------------------------------------------
    subroutine system_update_quantity(this, iq)
      import system_t
      import clock_t
      class(system_t),      intent(inout) :: this
      integer,              intent(in)    :: iq
    end subroutine system_update_quantity
 
    ! ---------------------------------------------------------
    subroutine system_update_exposed_quantity(partner, iq)
      import system_t
      import clock_t
      class(system_t), intent(inout) :: partner
      integer,         intent(in)    :: iq
    end subroutine system_update_exposed_quantity
 
    ! ---------------------------------------------------------
    subroutine system_restart_write_data(this)
      import system_t
      class(system_t), intent(inout) :: this
    end subroutine system_restart_write_data
 
    ! ---------------------------------------------------------
    ! this function returns true if restart data could be read
    logical function system_restart_read_data(this)
      import system_t
      class(system_t), intent(inout) :: this
    end function system_restart_read_data
    subroutine system_update_kinetic_energy(this)
      import system_t
      class(system_t),      intent(inout) :: this
    end subroutine system_update_kinetic_energy
 
  end interface
 
  type, extends(partner_list_t) :: system_list_t
    private
  contains
    procedure :: add => system_list_add_node 
    procedure :: contains => system_list_contains 
  end type system_list_t
 
  type, extends(linked_list_iterator_t) :: system_iterator_t
    private
  contains
    procedure :: get_next => system_iterator_get_next 
  end type system_iterator_t
 
contains
 
  ! ---------------------------------------------------------
  subroutine system_execute_algorithm(this)
    class(system_t),     intent(inout) :: this
 
    type(algorithmic_operation_t) :: operation
    logical :: all_updated, at_barrier, operation_done
    type(event_handle_t) :: debug_handle
    integer :: iq, iuq
    integer, allocatable :: updated_quantities(:)
 
    push_sub(system_execute_algorithm)
 
    at_barrier = .false.
 
    do while (.not. at_barrier)
 
      operation = this%algo%get_current_operation()
 
      debug_handle = multisystem_debug_write_event_in(this%namespace, event_function_call_t("dt_operation", operation),    &
        system_clock=this%clock, prop_clock=this%algo%clock)
 
      ! First try to execute the operation as a system specific operation
      operation_done = this%do_algorithmic_operation(operation, updated_quantities)
      if (allocated(updated_quantities)) then
        ! Update the quantities clock
        do iuq = 1, size(updated_quantities)
          iq = updated_quantities(iuq)
          if (.not. this%algo%inside_scf) then
            this%quantities(iq)%clock = this%quantities(iq)%clock + clock_tick
          else
            ! We set it to the propagation time to avoid double increment
            call this%quantities(iq)%clock%set_time(this%algo%clock)
          end if
 
          call multisystem_debug_write_marker(this%namespace, event_clock_update_t("quantity", quantity_label(iq), &
            this%quantities(iq)%clock, "set"))
        end do
      end if
 
      ! If not done, we try to execute it as an algorithm-specific operation.
      if (.not. operation_done) then
        operation_done = this%algo%do_operation(operation)
      else
        call this%algo%next()
      end if
 
      ! If still not done, the operation must be a generic operation
      if (.not. operation_done) then
 
        select case (operation%id)
        case (skip)
          ! Do nothing
          call this%algo%next()
 
        case (step_done)
          ! Increment the system clock by one time-step
          this%clock = this%clock + clock_tick
          call multisystem_debug_write_marker(this%namespace, event_clock_update_t("system",  "", this%clock, "tick"))
 
          ! Recompute the total energy
          call this%update_total_energy()
 
          ! Write output
          call this%output_write()
 
          ! Update elapsed time
          call this%algo%update_elapsed_time()
 
          ! Print information about the current iteration
          ! (NB: needs to be done after marking the propagation step as finished,
          ! so that the timings are correct)
          call this%iteration_info()
 
          call this%algo%next()
 
        case (rewind_algorithm)
          if (.not. this%arrived_at_any_barrier() .and. .not. this%algorithm_finished()) then
            ! Reset propagator for next step if not waiting at barrier and if
            ! the algorithm is not finished
            call this%algo%rewind()
          else
            at_barrier = .true.
          end if
 
        case (update_interactions)
          ! We increment by one algorithmic step
          this%algo%clock = this%algo%clock + clock_tick
          call multisystem_debug_write_marker(this%namespace, event_clock_update_t("propagator", "", this%algo%clock, "tick"))
 
          ! Try to update all the interactions
          all_updated = this%update_interactions()
 
          ! Move to next algorithm step if all interactions have been
          ! updated. Otherwise try again later.
          if (all_updated) then
            this%accumulated_loop_ticks = this%accumulated_loop_ticks + 1
            call this%algo%next()
          else
            this%algo%clock = this%algo%clock - clock_tick
            call multisystem_debug_write_marker(this%namespace, event_clock_update_t("propagator", "", this%algo%clock, "reverse"))
          end if
 
          ! Interactions are implicit barriers
          at_barrier = .true.
 
        case default
          message(1) = "Unsupported algorithmic operation."
          write(message(2), '(A,A,A)') trim(operation%id), ": ", trim(operation%label)
          call messages_fatal(2, namespace=this%namespace)
        end select
      end if
 
      call multisystem_debug_write_event_out(debug_handle, system_clock=this%clock, prop_clock=this%algo%clock)
    end do
 
    pop_sub(system_execute_algorithm)
  end subroutine system_execute_algorithm
 
  ! ---------------------------------------------------------
  subroutine system_reset_clocks(this, accumulated_ticks)
    class(system_t),      intent(inout) :: this
    integer,              intent(in)    :: accumulated_ticks
 
    integer :: iq
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
 
    character(len=MAX_INFO_LEN) :: extended_label
 
    push_sub(system_reset_clocks)
 
    ! Propagator clock
    this%algo%clock = this%algo%clock - accumulated_ticks*clock_tick
    call multisystem_debug_write_marker(this%namespace, event_clock_update_t("propagator", "", this%algo%clock, "reset"))
 
    ! Interaction clocks
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      interaction%clock = interaction%clock - accumulated_ticks*clock_tick
 
      select type (interaction)
      class is (interaction_with_partner_t)
        extended_label = trim(interaction%label)//"-"//trim(interaction%partner%namespace%get())
      class default
        extended_label = trim(interaction%label)
      end select
      call multisystem_debug_write_marker(this%namespace, event_clock_update_t( "interaction", extended_label, &
        interaction%clock, "reset"))
    end do
 
    ! Internal quantities clocks
    do iq = 1, max_quantities
      if (this%quantities(iq)%required) then
        this%quantities(iq)%clock = this%quantities(iq)%clock - accumulated_ticks*clock_tick
        call multisystem_debug_write_marker(this%namespace, event_clock_update_t("quantity", quantity_label(iq), &
          this%quantities(iq)%clock, "reset"))
      end if
    end do
 
    pop_sub(system_reset_clocks)
  end subroutine system_reset_clocks
 
  ! ---------------------------------------------------------
  logical function system_update_exposed_quantities(partner, requested_time, interaction) result(allowed_to_update)
    class(system_t),      intent(inout) :: partner
    type(clock_t),        intent(in)    :: requested_time
    class(interaction_t), intent(inout) :: interaction
 
    logical :: ahead_in_time, right_on_time, need_to_copy
    integer :: iq, q_id
 
    type(event_handle_t) :: debug_handle
 
    push_sub(system_update_exposed_quantities)
 
    debug_handle = multisystem_debug_write_event_in(system_namespace = partner%namespace, &
      event = event_function_call_t("system_update_exposed_quantities"), &
      partner_clock = partner%clock, &
      requested_clock = requested_time, &
      interaction_clock = interaction%clock)
 
    select type (interaction)
    class is (interaction_with_partner_t)
 
      if (partner%algo%inside_scf .or. partner%algo%clock + clock_tick < requested_time) then
        ! we are inside an SCF cycle and therefore are not allowed to expose any quantities.
        ! or we are too much behind the requested time
        allowed_to_update = .false.
      else
        allowed_to_update = .true.
        need_to_copy = .true.
        do iq = 1, size(interaction%partner_quantities)
          ! Get the requested quantity ID
          q_id = interaction%partner_quantities(iq)
 
          ! All needed quantities must have been marked as required. If not, then fix your code!
          assert(partner%quantities(q_id)%required)
 
          ! First update the exposed quantities that are updated on demand
          if (partner%quantities(q_id)%updated_on_demand) then
            if(partner%quantities(q_id)%clock < requested_time .and. &
              (partner%quantities(q_id)%clock + clock_tick <= requested_time .or. &
              (partner%interaction_timing == option__interactiontiming__timing_retarded .and. &
              partner%quantities(q_id)%clock + clock_tick > requested_time))) then
              ! We can update because the partner will reach this time in the next sub-timestep
              ! For retarded interactions, we need to allow the quantity to get ahead by one step
              call partner%update_exposed_quantity(q_id)
 
              partner%quantities(q_id)%clock = partner%quantities(q_id)%clock + clock_tick
              call multisystem_debug_write_marker(partner%namespace, event_clock_update_t("quantity", quantity_label(q_id), &
                partner%quantities(q_id)%clock, "set"))
            end if
          end if
 
          if (partner%quantities(q_id)%available_at_any_time) then
            ! We ignore the clock times when a quantity is available at any time
            ahead_in_time = .false.
            right_on_time = .true.
          else
            ! Compare the clock times
            ahead_in_time = partner%quantities(q_id)%clock > requested_time
            right_on_time = partner%quantities(q_id)%clock == requested_time
          end if
 
          select case (partner%interaction_timing)
          case (option__interactiontiming__timing_exact)
            ! only allow interaction at exactly the same time
            allowed_to_update = allowed_to_update .and. right_on_time
            need_to_copy = allowed_to_update
          case (option__interactiontiming__timing_retarded)
            ! allow retarded interaction
            allowed_to_update = allowed_to_update .and. &
              (right_on_time .or. ahead_in_time)
            need_to_copy = need_to_copy .and. .not. ahead_in_time
          case default
            call messages_not_implemented("Method for interaction quantity timing", namespace=partner%namespace)
          end select
 
        end do
 
        ! If the quantities have been updated, we copy them to the interaction
        if (allowed_to_update .and. need_to_copy) then
          select type (interaction)
          type is (ghost_interaction_t)
            ! Nothing to copy. We still need to check that we are at the right
            ! time for the update though!
          class default
            call partner%copy_quantities_to_interaction(interaction)
          end select
        end if
      end if
 
    class default
      message(1) = "A system can only expose quantities to an interaction as a partner."
      call messages_fatal(1, namespace=partner%namespace)
    end select
 
    call multisystem_debug_write_event_out(debug_handle, update=allowed_to_update, &
      partner_clock = partner%clock, &
      requested_clock = requested_time, &
      interaction_clock = interaction%clock)
 
    pop_sub(system_update_exposed_quantities)
  end function system_update_exposed_quantities
 
  ! ---------------------------------------------------------
  subroutine system_init_all_interactions(this)
    class(system_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
 
    push_sub(system_init_all_interactions)
 
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      select type (interaction)
      type is (ghost_interaction_t)
        ! Skip the ghost interactions
      class is (interaction_with_partner_t)
        call this%init_interaction(interaction)
        call interaction%partner%init_interaction_as_partner(interaction)
      class default
        call this%init_interaction(interaction)
      end select
    end do
 
    pop_sub(system_init_all_interactions)
  end subroutine system_init_all_interactions
 
  ! ---------------------------------------------------------
  logical function system_update_interactions(this) result(all_updated)
    class(system_t),      intent(inout) :: this
 
    logical :: none_updated
    integer :: iq, q_id, n_quantities
    class(interaction_t), pointer :: interaction
    type(interaction_iterator_t) :: iter
 
    push_sub(system_update_interactions)
 
    ! Some systems might need to perform some specific operations before the
    ! update. This should only be done if no interaction has been updated yet,
    ! so that it is only done once.
    none_updated = .true.
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      if (interaction%clock == this%algo%clock) then
        none_updated = .false.
        exit
      end if
    end do
    if (none_updated) then
      call this%update_interactions_start()
    end if
 
    !Loop over all interactions
    all_updated = .true.
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
 
      if (.not. interaction%clock == this%algo%clock) then
        if  (allocated(interaction%system_quantities)) then
          n_quantities = size(interaction%system_quantities)
        else
          n_quantities = 0
        end if
 
        ! Update the system quantities that will be needed for computing the interaction
        do iq = 1, n_quantities
          ! Get requested quantity ID
          q_id = interaction%system_quantities(iq)
 
          ! All needed quantities must have been marked as required. If not, then fix your code!
          assert(this%quantities(q_id)%required)
 
          if (.not. this%quantities(q_id)%updated_on_demand) then
            ! Quantity must be at the correct time, otherwise the algorithm and
            ! the interaction are incompatible
            if (.not. this%quantities(q_id)%clock == this%algo%clock .and. .not. this%quantities(q_id)%available_at_any_time) then
              write(message(1), '(5a)') "Interaction ", trim(interaction%label), " is incompatible with the selected algorithm."
              write(message(2), '(3a)') "The interaction requires the ", trim(quantity_label(q_id)), &
                " at a time it is not available."
              call messages_fatal(2, namespace=this%namespace)
            end if
 
            ! We do not need to update quantities that are not updated on
            ! demand, as the algorithm takes care of doing that, so move on to
            ! next quantity
            cycle
          end if
 
          if (.not. this%quantities(q_id)%clock == this%algo%clock) then
            ! The requested quantity is not at the requested time, so we try to update it
 
            ! Sanity check: it should never happen that the quantity is in advance
            ! with respect to the requested time.
            if (this%quantities(q_id)%clock > this%algo%clock) then
              message(1) = "The quantity clock is in advance compared to the requested time."
              call messages_fatal(1, namespace=this%namespace)
            end if
 
            call this%update_quantity(q_id)
            call multisystem_debug_write_marker(this%namespace, event_clock_update_t("quantity", quantity_label(q_id), &
              this%quantities(q_id)%clock, "set"))
          end if
 
        end do
 
        ! We can now try to update the interaction
        all_updated = interaction%update(this%algo%clock) .and. all_updated
      end if
    end do
 
    ! Some systems might need to perform some specific operations after all the
    ! interactions have been updated
    if (all_updated) then
      call this%update_interactions_finish()
    end if
 
    pop_sub(system_update_interactions)
  end function system_update_interactions
 
  ! ---------------------------------------------------------
  subroutine system_update_interactions_start(this)
    class(system_t), intent(inout) :: this
 
    push_sub(system_update_interactions_start)
 
    ! By default nothing is done just before updating the interactions. Child
    ! classes that wish to change this behaviour should override this method.
 
    pop_sub(system_update_interactions_start)
  end subroutine system_update_interactions_start
 
  ! ---------------------------------------------------------
  subroutine system_update_interactions_finish(this)
    class(system_t), intent(inout) :: this
 
    push_sub(system_update_interactions_finish)
 
    ! By default nothing is done just after updating the interactions. Child
    ! classes that wish to change this behaviour should override this method.
 
    pop_sub(system_update_interactions_finish)
  end subroutine system_update_interactions_finish
 
  ! ---------------------------------------------------------
  subroutine system_restart_write(this)
    class(system_t), intent(inout) :: this
 
    logical :: restart_write
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
    integer :: ii
 
    push_sub(system_restart_write)
 
    call parse_variable(this%namespace, 'RestartWrite', .true., restart_write)
 
    if (restart_write) then
      ! do some generic restart steps here
      ! write clock restart data
      call this%clock%restart_write('restart_clock_system', this%namespace)
      call this%algo%clock%restart_write('restart_clock_propagator', this%namespace)
      call iter%start(this%interactions)
      do while (iter%has_next())
        interaction => iter%get_next()
        call interaction%restart_write(this%namespace)
      end do
      do ii = 1, max_quantities
        if (this%quantities(ii)%required) then
          call this%quantities(ii)%clock%restart_write('restart_clock_quantity_'//trim(quantity_label(ii)), &
            this%namespace)
        end if
      end do
      ! the following call is delegated to the corresponding system
      call this%restart_write_data()
      message(1) = "Wrote restart data for system "//trim(this%namespace%get())
      call messages_info(1, namespace=this%namespace)
    end if
 
    pop_sub(system_restart_write)
  end subroutine system_restart_write
 
  ! ---------------------------------------------------------
  ! this function returns true if restart data could be read
  logical function system_restart_read(this)
    class(system_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
    integer :: ii
 
    push_sub(system_restart_read)
 
    ! do some generic restart steps here
    ! read clock data
    system_restart_read = this%clock%restart_read('restart_clock_system', this%namespace)
    system_restart_read = system_restart_read .and. &
      this%algo%clock%restart_read('restart_clock_propagator', this%namespace)
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      system_restart_read = system_restart_read .and. interaction%restart_read(this%namespace)
      ! reduce by one because of the first UPDATE_INTERACTIONS
      interaction%clock = interaction%clock - clock_tick
    end do
    do ii = 1, max_quantities
      if (this%quantities(ii)%required) then
        system_restart_read = system_restart_read .and. &
          this%quantities(ii)%clock%restart_read('restart_clock_quantity_'//trim(quantity_label(ii)), &
          this%namespace)
      end if
      if (this%quantities(ii)%updated_on_demand) then
        ! decrease clock by one for on-demand quantities to account for initial update_interactions step
        this%quantities(ii)%clock = this%quantities(ii)%clock - clock_tick
      end if
    end do
    ! the following call is delegated to the corresponding system
    system_restart_read = system_restart_read .and. this%restart_read_data()
 
    if (system_restart_read) then
      message(1) = "Successfully read restart data for system "//trim(this%namespace%get())
      call messages_info(1, namespace=this%namespace)
    end if
 
    pop_sub(system_restart_read)
  end function system_restart_read
 
  ! ---------------------------------------------------------
  subroutine system_output_start(this)
    class(system_t), intent(inout) :: this
 
    push_sub(system_output_start)
 
    ! By default nothing is done to regarding output. Child classes that wish
    ! to change this behaviour should override this method.
 
    pop_sub(system_output_start)
  end subroutine system_output_start
 
  ! ---------------------------------------------------------
  subroutine system_output_write(this)
    class(system_t), intent(inout) :: this
 
    push_sub(system_output_write)
 
    ! By default nothing is done to regarding output. Child classes that wish
    ! to change this behaviour should override this method.
 
    pop_sub(system_output_write)
  end subroutine system_output_write
 
  ! ---------------------------------------------------------
  subroutine system_output_finish(this)
    class(system_t), intent(inout) :: this
 
    push_sub(system_output_finish)
 
    ! By default nothing is done to regarding output. Child classes that wish
    ! to change this behaviour should override this method.
 
    pop_sub(system_output_finish)
  end subroutine system_output_finish
 
  ! ---------------------------------------------------------
  subroutine system_init_algorithm(this, factory)
    class(system_t),            intent(inout) :: this
    class(algorithm_factory_t), intent(in)    :: factory
 
    integer :: ii
 
    push_sub(system_init_algorithm)
 
    call messages_experimental('Multi-system framework')
 
    this%algo => factory%create(this)
 
    call this%init_clocks()
 
    !%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('InteractionTiming', this%interaction_timing, namespace=this%namespace)
 
    do ii = 1, number_barriers
      this%barrier(ii)%active = .false.
      this%barrier(ii)%target_time = m_zero
    end do
 
    pop_sub(system_init_algorithm)
  end subroutine system_init_algorithm
 
  ! ---------------------------------------------------------------------------------------
  recursive function system_algorithm_finished(this) result(finished)
    class(system_t),       intent(in) :: this
    logical :: finished
 
    finished = this%algo%finished()
 
  end function system_algorithm_finished
 
  ! ---------------------------------------------------------
  subroutine system_init_clocks(this)
    class(system_t),            intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
 
    push_sub(system_init_clocks)
 
    ! Initialize propagator clock
    this%algo%clock = clock_t(time_step=this%algo%dt/this%algo%algo_steps)
 
    ! Initialize system clock
    this%clock = clock_t(time_step=this%algo%dt)
 
    ! Interactions clocks
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      interaction%clock = this%algo%clock - clock_tick
    end do
 
    ! Required quantities clocks
    where (this%quantities%required)
      this%quantities%clock = this%algo%clock
    end where
    where (this%quantities%updated_on_demand)
      this%quantities%clock = this%algo%clock - clock_tick
    end where
 
    pop_sub(system_init_clocks)
  end subroutine system_init_clocks
 
  ! ---------------------------------------------------------
  subroutine system_propagation_start(this)
    class(system_t),      intent(inout) :: this
 
    logical :: all_updated
    type(event_handle_t) :: debug_handle
    integer, allocatable :: updated_quantities(:)
 
    push_sub(system_propagation_start)
 
    debug_handle = multisystem_debug_write_event_in(this%namespace, event_function_call_t("system_propagation_start"), &
      system_clock = this%clock, prop_clock = this%algo%clock)
 
    ! Update interactions at initial time
    all_updated = this%update_interactions()
    if (.not. all_updated) then
      message(1) = "Unable to update interactions when initializing the propagation."
      call messages_fatal(1, namespace=this%namespace)
    end if
 
    ! System-specific and propagator-specific initialization step
    if (this%algo%start_step%id /= skip) then
      if (.not. this%do_algorithmic_operation(this%algo%start_step, updated_quantities)) then
        message(1) = "Unsupported algorithmic operation."
        write(message(2), '(A,A,A)') trim(this%algo%start_step%id), ": ", trim(this%algo%start_step%label)
        call messages_fatal(2, namespace=this%namespace)
      else if (allocated(updated_quantities)) then
        message(1) = "Update of quantities not allowed in algorithmic operation."
        write(message(2), '(A,A,A)') trim(this%algo%final_step%id), ": ", trim(this%algo%final_step%label)
        call messages_fatal(2, namespace=this%namespace)
      end if
    end if
 
    ! Compute the total energy at the beginning of the simulation
    call this%update_total_energy()
 
    ! Start output
    call this%output_start()
 
    ! Write header for propagation log
    call messages_print_with_emphasis(msg="Multi-system propagation", namespace=this%namespace)
    write(message(1), '(a6,1x,a14,1x,a13,1x,a10,1x,a15)') 'Iter', 'Time', 'Energy', 'SC Steps', 'Elapsed Time'
    call messages_info(1, namespace=this%namespace)
    call messages_print_with_emphasis(namespace=this%namespace)
 
    ! Rewind propagator (will also set the initial time to compute the elapsed time)
    call this%algo%rewind()
 
    call multisystem_debug_write_event_out(debug_handle, system_clock = this%clock, prop_clock = this%algo%clock)
 
    pop_sub(system_propagation_start)
  end subroutine system_propagation_start
 
  ! ---------------------------------------------------------
  subroutine system_propagation_finish(this)
    class(system_t),      intent(inout) :: this
    type(event_handle_t) :: debug_handle
 
    integer, allocatable :: updated_quantities(:)
 
    push_sub(system_propagation_finish)
 
    debug_handle = multisystem_debug_write_event_in(this%namespace, event_function_call_t("system_propagation_finish"), &
      system_clock = this%clock, prop_clock = this%algo%clock)
 
    ! Finish output
    call this%output_finish()
 
    ! System-specific and propagator-specific finalization step
    if (this%algo%final_step%id /= skip) then
      if (.not.  this%do_algorithmic_operation(this%algo%final_step, updated_quantities)) then
        message(1) = "Unsupported algorithmic operation."
        write(message(2), '(A,A,A)') trim(this%algo%final_step%id), ": ", trim(this%algo%final_step%label)
        call messages_fatal(2, namespace=this%namespace)
      else if (allocated(updated_quantities)) then
        message(1) = "Update of quantities not allowed in algorithmic operation."
        write(message(2), '(A,A,A)') trim(this%algo%final_step%id), ": ", trim(this%algo%final_step%label)
        call messages_fatal(2, namespace=this%namespace)
      end if
    end if
 
    call multisystem_debug_write_event_out(debug_handle, system_clock = this%clock, prop_clock = this%algo%clock)
 
    pop_sub(system_propagation_finish)
  end subroutine system_propagation_finish
 
  ! ---------------------------------------------------------
  subroutine system_iteration_info(this)
    class(system_t), intent(in) :: this
 
    float :: energy
    character(len=40) :: fmt
 
    push_sub(system_iteration_info)
 
    energy = units_from_atomic(units_out%energy, this%total_energy)
    if (abs(energy) >= 1e5) then
      fmt = '(i7,1x,f14.6,1X,es13.6,1X,i9,1X,'
    else
      fmt = '(i7,1x,f14.6,1X,f13.6,1X,i9,1X,'
    end if
    if (this%algo%elapsed_time < 1e-3) then
      fmt = trim(fmt)//'es13.3)'
    else
      fmt = trim(fmt)//'f13.3)'
    end if
 
    write(message(1), fmt) this%clock%get_tick(), &
      units_from_atomic(units_out%time, this%clock%time()), energy, &
      0, this%algo%elapsed_time
    call messages_info(1, namespace=this%namespace)
 
    pop_sub(system_iteration_info)
  end subroutine system_iteration_info
 
  ! ---------------------------------------------------------
  logical function system_process_is_slave(this)
    class(system_t), intent(in) :: this
 
    push_sub(system_process_is_slave)
 
    ! By default an MPI process is not a slave
    system_process_is_slave = .false.
 
    pop_sub(system_process_is_slave)
  end function system_process_is_slave
 
  ! ---------------------------------------------------------
  subroutine system_end(this)
    class(system_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
 
    push_sub(system_end)
 
    ! No call to safe_deallocate macro here, as it gives an ICE with gfortran
    if (associated(this%algo)) then
      deallocate(this%algo)
    end if
 
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      safe_deallocate_p(interaction)
    end do
 
    pop_sub(system_end)
  end subroutine system_end
 
  ! ---------------------------------------------------------
  subroutine system_list_add_node(this, partner)
    class(system_list_t)         :: this
    class(interaction_partner_t), target :: partner
 
    push_sub(system_list_add_node)
 
    select type (partner)
    class is (system_t)
      call this%add_ptr(partner)
    class default
      assert(.false.)
    end select
 
    pop_sub(system_list_add_node)
  end subroutine system_list_add_node
 
  ! ---------------------------------------------------------
  recursive logical function system_list_contains(this, partner) result(contains)
    class(system_list_t)         :: this
    class(interaction_partner_t), target :: partner
 
    type(partner_iterator_t)  :: iterator
    class(interaction_partner_t),  pointer :: system
 
    push_sub(system_list_contains)
 
    contains = .false.
 
    select type (partner)
    class is (system_t)
 
      call iterator%start(this)
      do while (iterator%has_next() .and. .not. contains)
        system => iterator%get_next()
        contains = associated(system, partner)
      end do
 
    class default
      contains = .false.
    end select
 
    pop_sub(system_list_contains)
  end function system_list_contains
 
  ! ---------------------------------------------------------
  function system_iterator_get_next(this) result(system)
    class(system_iterator_t), intent(inout) :: this
    class(system_t),          pointer       :: system
 
    push_sub(system_iterator_get_next)
 
    select type (ptr => this%get_next_ptr())
    class is (system_t)
      system => ptr
    class default
      assert(.false.)
    end select
 
    pop_sub(system_iterator_get_next)
  end function system_iterator_get_next
 
  ! ---------------------------------------------------------
  subroutine system_init_parallelization(this, grp)
    class(system_t), intent(inout) :: this
    type(mpi_grp_t), intent(in)    :: grp
 
    push_sub(system_init_parallelization)
 
    call mpi_grp_copy(this%grp, grp)
    call messages_update_mpi_grp(this%namespace, grp)
 
    pop_sub(system_init_parallelization)
  end subroutine system_init_parallelization
 
 
 
  ! ---------------------------------------------------------
  subroutine system_start_barrier(this, target_time, barrier_index)
    class(system_t), intent(inout) :: this
    float,           intent(in)    :: target_time
    integer,         intent(in)    :: barrier_index
 
    push_sub(system_start_barrier)
 
    this%barrier(barrier_index)%active = .true.
    this%barrier(barrier_index)%target_time = target_time
 
    pop_sub(system_start_barrier)
  end subroutine system_start_barrier
 
  ! ---------------------------------------------------------
  subroutine system_end_barrier(this, barrier_index)
    class(system_t), intent(inout) :: this
    integer,         intent(in)    :: barrier_index
 
    push_sub(system_end_barrier)
 
    this%barrier(barrier_index)%active = .false.
    this%barrier(barrier_index)%target_time = m_zero
 
    pop_sub(system_end_barrier)
  end subroutine system_end_barrier
 
  ! ---------------------------------------------------------
  logical function system_arrived_at_barrier(this, barrier_index)
    class(system_t), intent(inout) :: this
    integer,         intent(in)    :: barrier_index
 
    type(clock_t) :: clock
 
    push_sub(system_arrived_at_barrier)
 
    system_arrived_at_barrier = .false.
    if (this%barrier(barrier_index)%active) then
      clock = this%clock + clock_tick
      if (clock%time() > this%barrier(barrier_index)%target_time) then
        system_arrived_at_barrier = .true.
      end if
    end if
 
    pop_sub(system_arrived_at_barrier)
  end function system_arrived_at_barrier
 
  ! ---------------------------------------------------------
  logical function system_arrived_at_any_barrier(this)
    class(system_t), intent(inout) :: this
 
    integer :: ii
 
    push_sub(system_arrived_at_any_barrier)
 
    system_arrived_at_any_barrier = .false.
    do ii = 1, number_barriers
      system_arrived_at_any_barrier = system_arrived_at_any_barrier &
        .or. this%arrived_at_barrier(ii)
    end do
 
    pop_sub(system_arrived_at_any_barrier)
  end function system_arrived_at_any_barrier
 
 
  ! ---------------------------------------------------------
  subroutine system_update_potential_energy(this)
    class(system_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
 
    push_sub(system_update_potential_energy)
 
    this%potential_energy = m_zero
 
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      if(.not. interaction%intra_interaction) then
        call interaction%calculate_energy()
        this%potential_energy = this%potential_energy + interaction%energy
      end if
    end do
 
    pop_sub(system_update_potential_energy)
  end subroutine system_update_potential_energy
 
  ! ---------------------------------------------------------
  subroutine system_update_internal_energy(this)
    class(system_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t), pointer :: interaction
 
    push_sub(system_update_internal_energy)
 
    this%internal_energy = m_zero
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      if(interaction%intra_interaction) then
        call interaction%calculate_energy()
        this%internal_energy = this%internal_energy + interaction%energy
      end if
    end do
 
    pop_sub(system_update_internal_energy)
  end subroutine system_update_internal_energy
 
  ! ---------------------------------------------------------
  subroutine system_update_total_energy(this)
    class(system_t), intent(inout) :: this
 
    push_sub(system_update_total_energy)
 
    call this%update_kinetic_energy()
    this%total_energy = this%kinetic_energy
 
    call this%update_potential_energy()
    this%total_energy = this%total_energy + this%potential_energy
 
    call this%update_internal_energy()
    this%total_energy = this%total_energy + this%internal_energy
 
    pop_sub(system_update_total_energy)
  end subroutine system_update_total_energy
 
end module system_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: