anchor.diff:     call mpi_comm_rank(mpi_comm_world, rank, err)
   
      !---------------------------------------------
  +   ! Compute number of domains
  +   !---------------------------------------------
  +   domain_size = ceiling(npoints * 1.0 / nprocs)
  +   N = ceiling(1._xp * npoints / nprocs)
  +
  +   !---------------------------------------------
      ! Read initial positions
      !---------------------------------------------
      open(newunit=un, file='initial_conditions.dat', status="old")
  -   call read_mpos(un, m, r)
  +   call read_mpos(un, rank*N + 1, (rank+1)*N + 1, m, r)
      close(un)
  +
  +   !---------------------------------------------
  +   ! Compute initial speeds and accelerations
  +   !---------------------------------------------
  +   call initial_speeds(r, v)
  +   istart = 1
  +   iend = npoints
  +   select case (flag_compute_force)
  +      case(0)
  +         call compute_force(m, r, istart, iend, a)
  +      case(1)
  +         call compute_force_omp(m, r, istart, iend, a)
  +      case(2)
  +         call compute_force_omp_nn_1(m, r, istart, iend, a)
  +      case default
  +         stop "Unknown value of flag_compute_force"
  +   end select
   
      !---------------------------------------------
      ! Print parameters
  
positive.code:             stop "Unknown value of flag_compute_mpi"
        end select
        !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ FIX BELOW
  
        !--------------------------------
        ! reduce accelerations
        !--------------------------------
  
negative.code: 
        call integrate(v, a, dt/2)  ! Compute v(t+dt/2)
        call integrate(r, v, dt)    ! Compute r(t+dt)
  
        select case(flag_compute_mpi)
           case(0)
  