anchor.diff:    integer :: error
   
     double precision :: sigma_N, sigma_C, max_cut
  -  double precision :: epsilon(3,2)
  +  double precision :: epsilon(3,2), shift
     double precision :: sigma(3,2), sigma_cut(3,2)
     double precision :: mass(2)
   
  
positive.code:     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
  
negative.code:      if (i .gt. 2000) then
          fixed = .false.
       end if
       
    end do setup
  
    call thermo_data%append(hfile, temperature, e1+e2+e_wall, kin_e, e1+e2+e_wall+kin_e, add=.false., force=.true.)
  