anchor.diff:            (solvent%Nmax + mass(1) + mass(2))
        call thermo_data%append(hfile, temperature, e1+e2+e_wall, kin_e, e1+e2+e_wall+kin_e, v_com)
   
  +     n_solvent = 0
  +     do k = 1, solvent%Nmax
  +        m = solvent%species(k)
  +        if (m <= 0) continue
  +        n_solvent(m) = n_solvent(m) + 1
  +     end do
  +     call n_solvent_el%append(n_solvent)
  +     call catalytic_change_el%append(catalytic_change)
  +     call bulk_change_el%append(bulk_change)
  +
        call dimer_io%position%append(colloids%pos)
        call dimer_io%velocity%append(colloids%vel)
        call dimer_io%image%append(colloids%image)
  
positive.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_loop*N_MD_steps, time=N_loop*N_MD_steps*dt)
    call h5gclose_f(fields_group, error)
  
    colloids% pos(3,:) = solvent_cells% edges(3)/2.d0
    if (order) then
  
negative.code:     integer :: m
      double precision :: x(3)
  
      catalytic_change = 0
      !$omp parallel do private(x, dist_to_C_sq, dist_to_N_sq) reduction(+:catalytic_change)
      do m = 1, solvent% Nmax
         if (solvent% flag(m) == 1) then
            x = rel_pos(colloids% pos(:,1), solvent% pos(:,m), solvent_cells% edges)
  