anchor.diff:    sigma(2,2) = 2*sigma_N
     sigma_cut = sigma*2**(1.d0/6.d0)
   
  -  call colloid_lj% init(epsilon, sigma, sigma_cut)
  +  call colloid_lj% init(epsilon(1:2,:), sigma(1:2,:), sigma_cut(1:2,:))
   
     epsilon(1,1) = 1.d0 
     epsilon(1,2) = 1.d0 
  
positive.code:      if (solvent% pos(2,m) < (L(2)/2.d0)) then
          solvent% species(m) = 1
       else
          solvent% species(m) = 2
       end if
    end do
  
  
negative.code:   sigma(2,:) = [sigma_C, sigma_N]
    sigma_cut = sigma*2**(1.d0/6.d0)
  
    call walls_colloid_lj% init(epsilon(1:2,:), sigma(1:2,:), sigma_cut(1:2,:))
  
    mass(1) = rho * sigma_C**3 * 4 * 3.14159265/3
    mass(2) = rho * sigma_N**3 * 4 * 3.14159265/3
    write(*,*) 'mass =', mass
  
    call solvent% init(N,3) !there will be 3 species of solvent particles
  
    call colloids% init(2,2, mass) !there will be 2 species of colloids
  
  