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 debug_oct_m
  use ghost_interaction_oct_m
  use global_oct_m
  use interaction_partner_oct_m
  use interaction_oct_m
  use iteration_counter_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 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_iteration_counters,                        &
    system_reset_iteration_counters,                       &
    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
    type(iteration_counter_t),    public :: iteration
    class(algorithm_t),  pointer, public :: algo => null()
 
    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_iteration_counters => system_reset_iteration_counters 
    procedure :: init_algorithm => system_init_algorithm 
    procedure :: algorithm_finished => system_algorithm_finished 
    procedure :: init_iteration_counters => system_init_iteration_counters 
    procedure :: init_all_interactions => system_init_all_interactions 
    procedure :: init_parallelization => system_init_parallelization 
    procedure :: update_couplings => system_update_couplings  
    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_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)
      use, intrinsic :: iso_fortran_env
      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_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_iteration=this%iteration, algo_iteration=this%algo%iteration)
 
      ! 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 iteration counters
        do iuq = 1, size(updated_quantities)
          iq = updated_quantities(iuq)
          call this%quantities(iq)%iteration%set(this%algo%iteration + 1)
          call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("quantity", quantity_label(iq), &
            this%quantities(iq)%iteration, "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 (iteration_done)
          ! Increment the system iteration by one
          this%iteration = this%iteration + 1
          call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("system",  "", this%iteration, "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_couplings)
          ! We increment by one algorithmic step
          this%algo%iteration = this%algo%iteration + 1
          call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("propagator", "", &
            this%algo%iteration, "tick"))
 
          ! Try to update all needed quantities from the interaction partners
          all_updated = this%update_couplings()
 
          ! Move to next algorithm step if all couplings have been
          ! updated. Otherwise try again later.
          if (all_updated) then
            call this%algo%next()
          else
            this%algo%iteration = this%algo%iteration - 1
            call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("propagator", "", &
              this%algo%iteration, "reverse"))
          end if
 
          ! Couplings are implicit barriers
          at_barrier = .true.
 
        case (update_interactions)
          ! Update all the interactions
          call this%update_interactions()
          call this%algo%next()
 
        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_iteration=this%iteration, algo_iteration=this%algo%iteration)
    end do
 
    pop_sub(system_execute_algorithm)
  end subroutine system_execute_algorithm
 
  ! ---------------------------------------------------------
  subroutine system_reset_iteration_counters(this, accumulated_iterations)
    class(system_t),      intent(inout) :: this
    integer,              intent(in)    :: accumulated_iterations
 
    integer :: iq
    type(interaction_iterator_t) :: iter
    class(interaction_t),              pointer :: interaction
 
    character(len=MAX_INFO_LEN) :: extended_label
 
    push_sub(system_reset_iteration_counters)
 
    ! Propagator counter
    this%algo%iteration = this%algo%iteration - accumulated_iterations
    call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("propagator", "", this%algo%iteration, "reset"))
 
    ! Interaction counters
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      interaction%iteration = interaction%iteration - accumulated_iterations
 
      extended_label = trim(interaction%label)//"-"//trim(interaction%partner%namespace%get())
      call multisystem_debug_write_marker(this%namespace, event_iteration_update_t( "interaction", extended_label, &
        interaction%iteration, "reset"))
    end do
 
    ! Internal quantities counters
    do iq = 1, max_quantities
      if (this%quantities(iq)%required) then
        this%quantities(iq)%iteration = this%quantities(iq)%iteration - accumulated_iterations
        call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("quantity", quantity_label(iq), &
          this%quantities(iq)%iteration, "reset"))
      end if
    end do
 
    pop_sub(system_reset_iteration_counters)
  end subroutine system_reset_iteration_counters
 
  ! ---------------------------------------------------------
  subroutine system_init_all_interactions(this)
    class(system_t), intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t),              pointer :: interaction
    integer :: timing
 
    push_sub(system_init_all_interactions)
 
    !%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 iteration or if retardation is allowed.
    !%Option timing_exact 1
    !% Only allow interactions at the exact iterations required by the algorithms behing executed
    !%Option timing_retarded 2
    !% Allow retarded interactions
    !%End
    call parse_variable(this%namespace, 'InteractionTiming', timing_exact, timing)
    if (.not. varinfo_valid_option('InteractionTiming', timing)) then
      call messages_input_error(this%namespace, 'InteractionTiming')
    end if
    call messages_print_var_option('InteractionTiming', timing, namespace=this%namespace)
 
    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 default
        call this%init_interaction(interaction)
        call interaction%partner%init_interaction_as_partner(interaction)
        interaction%timing = timing
 
        ! Sanity check for interaction timing and algorithm step
        if (interaction%timing == timing_exact) then
          select type (partner => interaction%partner)
          class is (system_t)
            if (this%algo%iteration%global_step() /= partner%algo%iteration%global_step() .and. &
              .not. all(partner%quantities(interaction%couplings_from_partner)%always_available)) then
              write(message(1), '(2a)') "InteractionTiming was set to exact timing, but systems ", &
                trim(this%namespace%get())
              write(message(2), '(3a)') "and ", trim(partner%namespace%get()), " have incompatible steps."
              call messages_fatal(2, namespace=this%namespace)
            end if
          end select
        end if
      end select
    end do
 
    pop_sub(system_init_all_interactions)
  end subroutine system_init_all_interactions
 
  ! ---------------------------------------------------------
  logical function system_update_couplings(this) result(all_updated)
    class(system_t),      intent(inout) :: this
 
    class(interaction_t),              pointer :: interaction
    type(interaction_iterator_t) :: iter
 
    push_sub(system_update_couplings)
 
    all_updated = .true.
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
 
      select type (partner => interaction%partner)
      class is (system_t)
        ! If the partner is a system, we can only update its couplings if it
        ! is not too much behind the requested iteration
        if (partner%algo%iteration + 1 >= this%algo%iteration) then
          call interaction%update_partner_couplings(this%algo%iteration)
        end if
 
      class default
        ! Partner that are not systems can have their couplings updated at any iteration
        call interaction%update_partner_couplings(this%algo%iteration)
      end select
 
      all_updated = all_updated .and. interaction%partner_couplings_up_to_date
    end do
 
    pop_sub(system_update_couplings)
  end function system_update_couplings
 
  ! ---------------------------------------------------------
  subroutine system_update_interactions(this)
    class(system_t),      intent(inout) :: this
 
    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.
    call this%update_interactions_start()
 
    !Loop over all interactions
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
 
      ! Update the system quantities that will be needed for computing the interaction
      if (allocated(interaction%system_quantities)) then
        n_quantities = size(interaction%system_quantities)
      else
        n_quantities = 0
      end if
      do iq = 1, n_quantities
        ! Get requested quantity ID
        q_id = interaction%system_quantities(iq)
 
        if (.not. this%quantities(q_id)%iteration == this%algo%iteration) then
          ! The requested quantity is not at the requested iteration, so we try to update it
 
          if (.not. this%quantities(q_id)%updated_on_demand) then
            ! Quantity must be at the correct iteration, otherwise the algorithm and
            ! the interaction are incompatible
            if (.not. this%quantities(q_id)%iteration == this%algo%iteration .and. &
              .not. this%quantities(q_id)%always_available) 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 iteration 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
 
          ! Sanity check: it should never happen that the quantity is in advance
          ! with respect to the requested iteration.
          if (this%quantities(q_id)%iteration > this%algo%iteration) then
            message(1) = "The quantity iteration is in advance compared to the requested iteration."
            call messages_fatal(1, namespace=this%namespace)
          end if
 
          call this%update_quantity(q_id)
          call multisystem_debug_write_marker(this%namespace, event_iteration_update_t("quantity", quantity_label(q_id), &
            this%quantities(q_id)%iteration, "set"))
        end if
 
      end do
 
      call interaction%update(this%algo%iteration)
    end do
 
    ! Some systems might need to perform some specific operations after all the
    ! interactions have been updated
    call this%update_interactions_finish()
 
    pop_sub(system_update_interactions)
  end subroutine 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 iteration counter restart data
      call this%iteration%restart_write('restart_iteration_system', this%namespace)
      call this%algo%iteration%restart_write('restart_iteration_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)%iteration%restart_write('restart_iteration_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 iteration data
    system_restart_read = this%iteration%restart_read('restart_iteration_system', this%namespace)
    system_restart_read = system_restart_read .and. &
      this%algo%iteration%restart_read('restart_iteration_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%iteration = interaction%iteration - 1
    end do
    do ii = 1, max_quantities
      if (this%quantities(ii)%required) then
        system_restart_read = system_restart_read .and. &
          this%quantities(ii)%iteration%restart_read('restart_iteration_quantity_'//trim(quantity_label(ii)), &
          this%namespace)
      end if
      if (this%quantities(ii)%updated_on_demand) then
        ! decrease iteration by one for on-demand quantities to account for initial update_interactions step
        this%quantities(ii)%iteration = this%quantities(ii)%iteration - 1
      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_iteration_counters()
 
    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_iteration_counters(this)
    class(system_t),            intent(inout) :: this
 
    type(interaction_iterator_t) :: iter
    class(interaction_t),              pointer :: interaction
 
    push_sub(system_init_iteration_counters)
 
    ! Initialize algorithm and system counters
    call this%algo%init_iteration_counters()
 
    ! Interactions counters
    call iter%start(this%interactions)
    do while (iter%has_next())
      interaction => iter%get_next()
      interaction%iteration = this%algo%iteration - 1
    end do
 
    ! Required quantities counters
    where (this%quantities%required)
      this%quantities%iteration = this%algo%iteration
    end where
    where (this%quantities%updated_on_demand)
      this%quantities%iteration = this%algo%iteration - 1
    end where
 
    pop_sub(system_init_iteration_counters)
  end subroutine system_init_iteration_counters
 
  ! ---------------------------------------------------------
  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_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    ! Update interactions at initial iteration
    all_updated = this%update_couplings()
    if (.not. all_updated) then
      message(1) = "Unable to update interactions when initializing the propagation."
      call messages_fatal(1, namespace=this%namespace)
    end if
    call this%update_interactions()
 
    ! System-specific and propagator-specific initialization step
    if (this%algo%start_operation%id /= skip) then
      if (.not. this%do_algorithmic_operation(this%algo%start_operation, updated_quantities)) then
        message(1) = "Unsupported algorithmic operation."
        write(message(2), '(A,A,A)') trim(this%algo%start_operation%id), ": ", trim(this%algo%start_operation%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%start_operation%id), ": ", trim(this%algo%start_operation%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_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    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_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    ! Finish output
    call this%output_finish()
 
    ! System-specific and propagator-specific finalization step
    if (this%algo%final_operation%id /= skip) then
      if (.not.  this%do_algorithmic_operation(this%algo%final_operation, updated_quantities)) then
        message(1) = "Unsupported algorithmic operation."
        write(message(2), '(A,A,A)') trim(this%algo%final_operation%id), ": ", trim(this%algo%final_operation%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_operation%id), ": ", trim(this%algo%final_operation%label)
        call messages_fatal(2, namespace=this%namespace)
      end if
    end if
 
    call multisystem_debug_write_event_out(debug_handle, system_iteration = this%iteration, algo_iteration = this%algo%iteration)
 
    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%iteration%counter(), &
      units_from_atomic(units_out%time, this%iteration%value()), 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()
      if (associated(interaction)) then
        deallocate(interaction)
      end if
    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(iteration_counter_t) :: iteration
 
    push_sub(system_arrived_at_barrier)
 
    system_arrived_at_barrier = .false.
    if (this%barrier(barrier_index)%active) then
      iteration = this%iteration + 1
      if (iteration%value() > 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: