anchor.diff:  
         integer :: err = 0
         integer :: stat(MPI_STATUS_SIZE)
  -      integer :: i                          ! Counter for elements
  +      integer :: i, np_i               ! Counters for elements
   
         integer :: istart, iend
   
  
positive.code: 
              do i = 0, nprocs / 2 - 1
                 ! Get data from nprocs-i in i
                 ! Note: use nprocs-i, because nprocs start at 0
  
                 if (rank == nprocs-i-1) then
                    call mpi_send(r, 3*N, MPI_REAL_XP, i, 0, MPI_COMM_WORLD, err)
                 else if (rank == i) then
                    call mpi_recv(r_right, 3*N, MPI_REAL_XP, nprocs-i-1, 0, MPI_COMM_WORLD, stat, err)
                 end if
  
                 if (rank == i) then
                    r_i = r
                 else if (rank == nprocs - i - 1) then
                    r_np_i = r
                 end if
  
                 ! Broadcast i-th data
                 call mpi_bcast(r_i, 3*N, MPI_REAL_XP, i, MPI_COMM_WORLD, err)
  
                 ! Broadcast n-i-th data
                 call mpi_bcast(r_np_i, 3*N, MPI_REAL_XP, nprocs-i-1, MPI_COMM_WORLD, err)
  
                 ! compute own interactions
                 a_comm_i = 0._xp
                 a_comm_np_i = 0._xp
                 if (rank == i) then
                    a_right = 0._xp
                    call compute_force_diag(m, r, 1, N, N, a)
                    call compute_force_diag(m, r_right, 1, N, N, a_right)
                    a_comm_i = a
  
                    ! compute interaction on right side
                 else if (rank < i) then
                    call compute_force(m, r_np_i, 1, N, r_right, 1, N, a_comm_np_i)
                    a_right = a_right - a_comm_np_i
  
                    ! compute interaction on left side
                 else
                    call compute_force(m, r_i, 1, N, r, 1, N, a_comm_i)
                    a = a - a_comm_i ! FIXME: False with different masses
  
                    if (rank == nprocs - i - 1) then
                       a_comm_np_i = a
                    end if
  
                 end if
  
                 ! Receive interaction in i
                 call mpi_reduce(a_comm_i,    a, 3*N, MPI_REAL_XP, MPI_SUM,          i, MPI_COMM_WORLD, err)
  
                 ! Receive interaction in np-i-1
                 call mpi_reduce(a_comm_np_i, a, 3*N, MPI_REAL_XP, MPI_SUM, nprocs-i-1, MPI_COMM_WORLD, err)
  
              end do
  
  
negative.code:             ! Gather missing data
              !--------------------------------
              do i = 0, nprocs / 2 - 1
                 np_i = nprocs - i - 1
                 if (rank == i) then
                    call mpi_send(a_right,     3*N, MPI_REAL_XP, np_i, 0, MPI_COMM_WORLD, err)
                 else if (rank == np_i) then
                    call mpi_recv(a_comm_np_i, 3*N, MPI_REAL_XP,    i, 0, MPI_COMM_WORLD, stat, err)
                    a = a + a_comm_np_i
                 end if
              end do
  
              deallocate(a_comm_i)
              deallocate(a_comm_np_i)
              deallocate(a_right)
              deallocate(r_i)
              deallocate(r_np_i)
              deallocate(r_right)
  
           case default
              stop "Unknown value of flag_mpi"
        end select
  