anchor.diff:  
   contains
   
  -  subroutine flag_particles_nl
  +  subroutine flag_particles
       double precision :: dist_to_C_sq
  -    integer :: i, r, s
  +    integer :: i, r, s, thread_id
       double precision :: x(3)
   
  +    !$omp parallel private(thread_id)
  +    thread_id = omp_get_thread_num() + 1
  +    !$omp do private(i, s, r, x, dist_to_C_sq)
       do i = 1, colloids%Nmax
          if (colloids%species(i)==1) then
             do s = 1,neigh% n(i)
  
positive.code:            colloids% vel(:,k) = colloids% vel(:,k) + &
               dt * ( colloids% force(:,k) + colloids% force_old(:,k) ) / (2 * colloids%mass(colloids%species(k)))
          end do
  
          call rattle_body_pos(colloids, links, links_d, dt, solvent_cells% edges, 1d-9)
          call varia%tac()
  
          call time_flag%tic()
          call flag_particles_nl
          call time_flag%tac()
          call time_change%tic()
          call change_species
  
negative.code:   call h5md_write_attribute(radial_hist_el%id, 'dx', radial_hist%dx)
    call h5gclose_f(fields_group, error)
  
    call h5gcreate_f(hfile%id, 'connectivity', connectivity_group, error)
    call h5md_write_dataset(connectivity_group, 'janus_links', links-1)
    call h5dopen_f(connectivity_group, 'janus_links', tmp_id, error)
    call h5md_write_attribute(tmp_id, 'particles_group', 'janus')
    call h5dclose_f(tmp_id, error)
    call h5gclose_f(connectivity_group, error)
  
    call solvent% random_placement(solvent_cells% edges, colloids, solvent_colloid_lj)
  
    call solvent% sort(solvent_cells)
  