anchor.diff:    type(lj_params_t) :: solvent_colloid_lj
     type(lj_params_t) :: colloid_lj
   
  -  integer, parameter :: rho = 10
  +  integer, parameter :: rho = 9
     integer :: N
     integer :: error
     double precision, allocatable :: pos_rattle(:,:)
  
positive.code:   call solvent_cells%init(L, 1.d0)
    colloids% pos(:,1) = solvent_cells% edges/2.0
    colloids% pos(:,2) = solvent_cells% edges/2.0 
    colloids% pos(1,2) = colloids% pos(1,2) + 2.0*sigma + 0.5d0
    
    call solvent% random_placement(solvent_cells% edges, colloids, solvent_colloid_lj)
  
    call solvent% sort(solvent_cells)
  
    call neigh% init(colloids% Nmax, 300)
    call neigh% make_stencil(solvent_cells, 1.5d0)
  
    call neigh% update_list(colloids, solvent, 1.5d0, solvent_cells)
  
    tau = 0.1d0 !was 0.1d0 but we changed it to Kaprals' value 
    N_MD_steps = 100
  
negative.code:      if (solvent% flag(m) == 1) then
         dist_to_C_sq = dot_product(colloids% pos(:,1) - solvent% pos(:,m),colloids% pos(:,1) - solvent% pos(:,m))
         dist_to_N_sq = dot_product(colloids% pos(:,2) - solvent% pos(:,m),colloids% pos(:,2) - solvent% pos(:,m))
         if ((dist_to_C_sq > (sigma_cut)**2) .and. (dist_to_N_sq > (sigma_cut)**2)) then
           solvent% species(m) = 2
           solvent% flag(m) = 0
         end if
  