anchor.diff:      colloid_pos(:,2) = 0
       colloid_pos(3,2) = d + colloids% pos(3,1)
     end subroutine concentration_field
  -  
  -  function colloid_wall_interaction(colloids, walls_colloid_lj,edges) result(e)
  -     type(particle_system_t), intent(inout) :: colloids
  -     type(lj_params_t), intent(in) :: walls_colloid_lj
  -     double precision, intent(in) :: edges(3)
  -     double precision :: e
  -
  -     integer :: k
  -     integer :: s
  -     double precision :: r_sq1, d1(3), r_sq2, d2(3),f1(3),f2(3)
  -
  -     f1= 0.d0
  -     f2 = 0.d0
  -     e = 0.d0
  -     do k = 1, colloids% Nmax
  -        s = colloids% species(k)
  -        d1 = rel_pos(colloids% pos(:,k), [edges(1)/2.d0, edges(2)/2.d0,0.d0],edges)
  -        d2 = rel_pos(colloids% pos(:,k), [edges(1)/2.d0, edges(2)/2.d0,edges(3)],edges)
  -        r_sq1 = dot_product(d1,d1)
  -        r_sq2 = dot_product(d2,d2)
  -        if (r_sq1 <= walls_colloid_lj% cut_sq(1,s)) then
  -           f1 = lj_force_9_6(d1, r_sq1, walls_colloid_lj% epsilon(1,s), walls_colloid_lj% sigma(1,s))
  -           e = e + lj_energy_9_6(r_sq1, walls_colloid_lj% epsilon(1,s), walls_colloid_lj% sigma(1,s))
  -        end if
  -        if (r_sq2 <= walls_colloid_lj% cut_sq(1,s)) then
  -           f2 = lj_force_9_6(d2, r_sq2, walls_colloid_lj% epsilon(1,s), walls_colloid_lj% sigma(1,s))
  -           e = e + lj_energy_9_6(r_sq2, walls_colloid_lj% epsilon(1,s), walls_colloid_lj% sigma(1,s))
  -        end if
  -        colloids% force(:,k) = colloids% force(:,k) - f1 - f2
  -     end do
  -
  -  end function colloid_wall_interaction 
   
     subroutine mpcd_stream_zwall_light(particles, cells, dt,g)
       type(particle_system_t), intent(inout) :: particles
  
positive.code: 
    e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
    e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
    !e_wall = colloid_wall_interaction(colloids, walls_colloid_lj,solvent_cells% edges)
    solvent% force_old = solvent% force
    colloids% force_old = colloids% force
  
  
negative.code:   integer :: error
  
    double precision :: sigma_N, sigma_C, max_cut
    double precision :: epsilon(3,2), shift
    double precision :: sigma(3,2), sigma_cut(3,2)
    double precision :: mass(2)
  
  