anchor.diff:        if (mod(iter, dump_freq) == 0) then
            print*, 'Dump!'
            call compute_energy(m, r, v, Ec, Ep, E) ! Compute Ec, Ep, E at t+dt
  -         call write_dump(11, 12, iter, Ec, Ep, E, t, r, v)
  +         call write_dump(una, un, iter, Ec, Ep, E, t, r, v)
         end if
   
      end do
   
  -   close(unit=11)
  -   close(unit=12)
  -
      !---------------------------------------------
  -   ! Deallocations
  +   ! Close files
      !---------------------------------------------
  -   deallocate(r)
  -   deallocate(v)
  -   deallocate(a)
  +   close(un)
  +   close(una)
   
      !---------------------------------------------
      ! Elapsed time
  
positive.code:    call compute_initial_variables()
  
     !---------------------------------------------
     ! Compute time step
     !---------------------------------------------
     ! the timestep here is experimental and corresponds to 1% of the
     ! dynamic time
     dt = 1.e-2_xp
     print *, dt
  
     !---------------------------------------------
     ! Loop over time
     !---------------------------------------------
     open(unit=11, file='output.dat')
     open(unit=12, file='output_int.dat')
  
     call write_dump_headers(12)
  
     do while (iter < maxiter)
  
        call integrate(v, a, dt/2)  ! Compute v(t+dt/2)
  
negative.code:    !---------------------------------------------
     call cpu_time(start_time)
  
     real(kind = xp) :: m(npoints)        ! Masses of the particles
     real(kind = xp) :: r(npoints,3)      ! Positions of the particles (3-dim vectors)
     real(kind = xp) :: v(npoints,3)      ! Speeds of the particles (3-dim vectors)
     real(kind = xp) :: a(npoints,3)      ! Acceleration of the particles (3-dim vectors)
     real(kind = xp) :: t = 0._xp         ! Total time elapsed in the simulation
     real(kind = xp) :: dt                ! Timestep
     real(kind = xp) :: Ec                ! Total kinetic energy
     real(kind = xp) :: Ep                ! Total potential energy
     real(kind = xp) :: E                 ! Total energy
     integer         :: iter      = 0     ! Number of iterations ran
     integer         :: dump_freq = 10    ! Frequency at which the system is sampled
     integer         :: maxiter   = 10000 ! Maximum number of iterations
  
     !---------------------------------------------
     ! Read initial positions
     !---------------------------------------------
     open(newunit=un, file='initial_conditions.dat', status="old")
     call read_mpos(un, m, r)
     close(un)
  
     !---------------------------------------------
     ! Compute initial speeds, accelerations, variables
  