anchor.diff:    wall_v = 0
     wall_t = [T, T]
     
  -  tau =0.1d0 !PTread_d(config, 'tau')
  +  tau =0.2d0 !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')
  +  N_loop = 20000 !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)
  +  epsilon(:,1) = 1 !PTread_dvec(config, 'epsilon_C', 2)
  +  epsilon(:,2) = [1.d0, 0.1d0, 1.d0] !PTread_dvec(config, 'epsilon_N', 2)
   
     sigma(:,1) = sigma_C
     sigma(:,2) = sigma_N
  
positive.code:   sigma(2,:) = [sigma_C, sigma_N]
    sigma_cut = sigma*2**(1.d0/6.d0)
  
    call walls_colloid_lj% init(epsilon, sigma, sigma_cut)
  
    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,2) !there will be 2 species of solvent particles
  
    call colloids% init(2,2, mass) !there will be 2 species of colloids
  
  
negative.code:   sigma(2,2) = 2*sigma_N
    sigma_cut = sigma*2**(1.d0/6.d0)
  
    call colloid_lj% init(epsilon(1:2,:), sigma(1:2,:), sigma_cut(1:2,:))
  
    epsilon(1,1) = 1.d0 
    epsilon(1,2) = 1.d0 
  