anchor.diff:  
     epsilon = 1.d0
     sigma(1,:) = [sigma_C, sigma_N]
  -  sigma(2,:) = [sigma_C, sigma_N]
     sigma_cut = sigma*3**(1.d0/6.d0)
  -  shift = maxval(colloid_lj% cut) +0.25
  -  call walls_colloid_lj% init(epsilon(1:2,:), sigma(1:2,:), sigma_cut(1:2,:), shift)
  +  shift = max(sigma_C, sigma_N)*2**(1./6.) + 0.25
  +  call walls_colloid_lj% init(epsilon(1:1,:), sigma(1:1,:), sigma_cut(1:1,:), shift)
  +  write(*,*) epsilon(1:2,:), sigma(1:2,:), sigma_cut(1:2,:), shift
   
     mass(1) = rho * sigma_C**3 * 4 * 3.14159265/3
     mass(2) = rho * sigma_N**3 * 4 * 3.14159265/3
  
positive.code:         call vx% reset()
       end if
  
       call solvent% sort(solvent_cells)
       call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells)
  
  
negative.code:   double precision :: colloid_pos(3,2)
    type(h5md_file_t) :: hfile
    type(h5md_element_t) :: dummy_element
    integer(HID_T) :: fields_group
    type(h5md_element_t) :: rho_xy_el
    type(thermo_t) :: thermo_data
    type(particle_system_io_t) :: dimer_io
    type(particle_system_io_t) :: solvent_io
  