anchor.diff:                 then
                solvent% species(m) = 2
                solvent% flag(m) = 0
  +             catalytic_change(1) = catalytic_change(1) - 1
  +             catalytic_change(2) = catalytic_change(2) + 1
             end if
          end if
       end do
  
positive.code:      double precision, intent(in) :: edges(3)
       integer, intent(in) :: bufferlength
    
       integer :: k  
  
       do k = 1, particles% Nmax
          if (particles% pos(1,k) < bufferlength) then
             if (particles% pos(2,k) < edges(2)/2.d0) then
                particles% species(k) = 1
             else
                particles% species(k) = 3
             end if  
          end if 
       end do
  
negative.code:   allocate(rho_xy(N_species, L(2), L(1)))
    call h5gcreate_f(hfile%id, 'fields', fields_group, error)
    call rho_xy_el%create_time(fields_group, 'rho_xy', rho_xy, ior(H5MD_LINEAR,H5MD_STORE_TIME), &
         step=N_MD_steps, time=N_MD_steps*dt)
    call h5gclose_f(fields_group, error)
  
    call n_solvent_el%create_time(hfile%observables, 'n_solvent', &
         n_solvent, ior(H5MD_LINEAR,H5MD_STORE_TIME), step=N_MD_steps, &
         time=N_MD_steps*dt)
    call catalytic_change_el%create_time(hfile%observables, 'catalytic_change', &
         catalytic_change, ior(H5MD_LINEAR,H5MD_STORE_TIME), step=N_MD_steps, &
         time=N_MD_steps*dt)
    call bulk_change_el%create_time(hfile%observables, 'bulk_change', &
         bulk_change, ior(H5MD_LINEAR,H5MD_STORE_TIME), step=N_MD_steps, &
         time=N_MD_steps*dt)
  
    colloids% pos(3,:) = solvent_cells% edges(3)/2.d0
    if (order) then
  