anchor.diff:                   colloids% vel
   
        call solvent% sort(solvent_cells)
  -     call neigh% update_list(colloids, solvent, 1.5d0, solvent_cells)
  +     call neigh% update_list(colloids, solvent, 1.5d0*sigma, solvent_cells)
   
        call simple_mpcd_step(solvent, solvent_cells, mt)
   
  
positive.code:   type(lj_params_t) :: solvent_colloid_lj
    type(lj_params_t) :: colloid_lj
  
    integer, parameter :: rho = 10
    integer :: N
    integer :: error
    double precision, allocatable :: pos_rattle(:,:)
  
negative.code: 
    call h5open_f(error)
  
    L = [32, 32, 32]
    N = rho *L(1)*L(2)*L(3)
    
    allocate(epsilon(2,2))
  
    epsilon = reshape((/ 1.d0, 1.d0, 1.d0, 0.5d0 /),shape(epsilon))
    sigma = 2.d0
    sigma_cut = sigma*2**(1.d0/6.d0)
  
    call solvent_colloid_lj% init( epsilon, &
         reshape( [ sigma, sigma, sigma, sigma ], [2, 2] ), reshape( [ sigma_cut,sigma_cut,sigma_cut,sigma_cut ], [2, 2] ) )
    
    epsilon1 = 1.d0
    sigma = 2.d0
    sigma_cut = sigma*2**(1.d0/6.d0)
  
    call colloid_lj% init( reshape( [ epsilon1,epsilon1,epsilon1,epsilon1 ], [2, 2] ), &
  