Fix OEP full - spin polarized case

Description

The OEP full was not taking the right occupations for the spin-polarized case.

News snippet

The OEP full was not taking the right occupations for the spin-polarized case.

Checklist

• I have checked that my code follows the Octopus coding standards
• I have added tests for all the new features added in this request.

src/multisystem/multisystem.F90

@@ -54,6 +54,7 @@ module multisystem_oct_m
     procedure :: dt_operation =>  multisystem_dt_operation
     procedure :: init_parallelization => multisystem_init_parallelization
     procedure :: smallest_algo_dt => multisystem_smallest_algo_dt
+    procedure :: largest_dt => multisystem_largest_dt 
     procedure :: next_time_on_largest_dt => multisystem_next_time_on_largest_dt
     procedure :: reset_clocks => multisystem_reset_clocks
     procedure :: init_propagator => multisystem_init_propagator
@@ -216,6 +217,33 @@ contains
     POP_SUB(multisystem_smallest_algo_dt)
   end function multisystem_smallest_algo_dt
 
+  ! ---------------------------------------------------------------------------------------
+  recursive function multisystem_largest_dt(this) result(largest_dt)
+    class(multisystem_t), intent(inout) :: this
+    FLOAT :: largest_dt
+
+    type(system_iterator_t) :: iter
+    class(system_t), pointer :: system
+
+    PUSH_SUB(multisystem_largest_dt)
+
+    largest_dt = M_ZERO
+    call iter%start(this%list)
+    do while (iter%has_next())
+      system => iter%get_next()
+      select type (system)
+      class is (multisystem_t)
+        largest_dt = max(largest_dt, system%largest_dt())
+      class default
+        largest_dt = max(largest_dt, system%prop%dt)
+      end select
+    end do
+
+    POP_SUB(multisystem_largest_dt)
+  end function multisystem_largest_dt
+
+
+
   ! ---------------------------------------------------------------------------------------
   recursive function multisystem_next_time_on_largest_dt(this) result(next_time_on_largest_dt)
     class(multisystem_t), intent(inout) :: this
@@ -312,6 +340,8 @@ contains
     class(system_t), pointer :: system
     type(interaction_iterator_t) :: inter_iter
     class(interaction_t), pointer :: interaction
+    FLOAT :: smallest_algo_dt, largest_dt
+    integer :: nsteps
 
     PUSH_SUB(multisystem_init_propagator)
 
@@ -322,13 +352,17 @@ contains
       call system%init_propagator()
     end do
 
+    smallest_algo_dt = this%smallest_algo_dt()
+    largest_dt = this%largest_dt()
+    nsteps = int(largest_dt/smallest_algo_dt)
+
     ! 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())
+    ! as we use the largest dt of the subsystems.
+    this%prop => propagator_static_t(largest_dt, nsteps)
     this%interaction_timing = OPTION__INTERACTIONTIMING__TIMING_EXACT
     call this%prop%rewind()