anchor.diff:       if (solvent% pos(2,m) < (L(2)/2.d0)) then
           solvent% species(m) = 1
        else
  -        solvent% species(m) = 2
  +        solvent% species(m) = 3
        end if
     end do
   
  
positive.code:   wall_v = 0
    wall_t = [T, T]
    
    tau =0.1d0 !PTread_d(config, 'tau')
    N_MD_steps = 100 !PTread_i(config, 'N_MD')
    dt = tau / N_MD_steps
    N_loop = 10 !PTread_i(config, 'N_loop')
  
    
  
    sigma_C = 2.d0 !PTread_d(config, 'sigma_C')
    sigma_N = 2.d0 !PTread_d(config, 'sigma_N')
    
    epsilon(:,1) = [1.d0, 0.1d0] !PTread_dvec(config, 'epsilon_C', 2)
    epsilon(:,2) = [1.d0, 1.d0] !PTread_dvec(config, 'epsilon_N', 2)
  
    sigma(:,1) = sigma_C
    sigma(:,2) = sigma_N
  
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
  
  