anchor.diff:    sigma(2,2) = 2*sigma_N
     sigma_cut = sigma*2**(1.d0/6.d0)
   
  -  
  -
     call colloid_lj% init(epsilon, sigma, sigma_cut)
   
     mass(1) = rho * sigma_C**3 * 4 * 3.14159265/3
  
positive.code: 
          do k=1, colloids% Nmax
             colloids% vel(:,k) = colloids% vel(:,k) + &
               dt * ( colloids% force(:,k) + colloids% force_old(:,k) ) / (2 * mass(colloids%species(k)))
          end do
  
          call rattle_dimer_vel(colloids, d, dt, solvent_cells% edges)
  
negative.code: 
    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
  
    call PTkill(config)
    
  