anchor.diff:    call threefry_rng_init(state, args%seed)
   
     call h5open_f(error)
  -
  -  L = PTread_ivec(config, 'L', 3)
  -  rho = PTread_i(config, 'rho')
  -  T_init = PTread_d(config, 'T')
  +  call hfile%create(args%output_file, 'RMPCDMD::n_colloids_pbc', &
  +       RMPCDMD_REVISION, 'Pierre de Buyl')
  +  call h5gcreate_f(hfile%id, 'parameters', params_group, error)
  +  call hdf5_util_write_dataset(params_group, 'seed', args%seed)
  +
  +  L = PTread_ivec(config, 'L', 3, loc=params_group)
  +  rho = PTread_i(config, 'rho', loc=params_group)
  +  T_init = PTread_d(config, 'T', loc=params_group)
     T = T_init
  -  T_final = PTread_d(config, 'T_final')
  -  tau = PTread_d(config, 'tau')
  -  N_MD_steps = PTread_i(config, 'N_MD')
  +  T_final = PTread_d(config, 'T_final', loc=params_group)
  +  tau = PTread_d(config, 'tau', loc=params_group)
  +  N_MD_steps = PTread_i(config, 'N_MD', loc=params_group)
  +  colloid_sampling = PTread_i(config, 'colloid_sampling', loc=params_group)
     dt = tau / N_MD_steps
  -  N_loop = PTread_i(config, 'N_loop')
  -
  -  N_colloids = PTread_i(config, 'N_colloids')
  -  epsilon = PTread_d(config, 'epsilon')
  -  sigma = PTread_d(config, 'sigma')
  +  N_loop = PTread_i(config, 'N_loop', loc=params_group)
  +
  +  N_colloids = PTread_i(config, 'N_colloids', loc=params_group)
  +  epsilon = PTread_d(config, 'epsilon', loc=params_group)
  +  sigma = PTread_d(config, 'sigma', loc=params_group)
     sigma_cut = sigma*2.d0**(1.d0/6.d0)
  +
  +  call thermo_data%init(hfile, n_buffer=50, step=N_MD_steps, time=N_MD_steps*dt)
   
     mass = rho * sigma**3 * 4 * 3.141/3
     call colloids%init(N_colloids, 1, [mass])
     colloids%species = 1
   
  +  colloids_io%force_info%store = .false.
  +  colloids_io%id_info%store = .false.
  +  colloids_io%velocity_info%store = .false.
  +  colloids_io%position_info%store = .true.
  +  colloids_io%position_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  +  colloids_io%position_info%step = N_MD_steps
  +  colloids_io%position_info%time = N_MD_steps*dt
  +  colloids_io%image_info%store = .true.
  +  colloids_io%image_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  +  colloids_io%image_info%step = N_MD_steps
  +  colloids_io%image_info%time = N_MD_steps*dt
  +  colloids_io%species_info%store = .true.
  +  colloids_io%species_info%mode = H5MD_FIXED
  +  call colloids_io%init(hfile, 'colloids', colloids)
  +
  +  call h5gcreate_f(colloids_io%group, 'box', box_group, error)
  +  call h5md_write_attribute(box_group, 'dimension', 3)
  +  call dummy_element%create_fixed(box_group, 'edges', solvent_cells%edges)
  +  call h5gclose_f(box_group, error)
  +
     N = rho*L(1)*L(2)*L(3) - int(rho*4*3.142/3 * sigma**3*colloids%Nmax)
  -  write(*,*) N, 'colloids'
  +
  +  write(*,*) N, 'solvent particles'
   
     if (N <= 0) error stop 'Not enough volume available for solvent'
   
     call solvent_colloid_lj% init( reshape( [ epsilon ], [1, 1] ), &
          reshape( [ sigma ], [1, 1] ), reshape( [ sigma_cut ], [1, 1] ) )
   
  -  epsilon = PTread_d(config, 'epsilon_colloids')
  +  epsilon = PTread_d(config, 'epsilon_colloids', loc=params_group)
     sigma_cut = 3*sigma
  +
  +  call h5gclose_f(params_group, error)
  +  call PTkill(config)
   
     call colloid_lj% init( reshape( [ epsilon ], [1, 1] ), &
          reshape( [ 2*sigma ], [1, 1] ), reshape( [ 2*sigma_cut ], [1, 1] ) )
  
positive.code:   double precision :: so_max, co_max
  
    double precision :: e1, e2
    double precision :: tau, dt, T, T_init, T_final
    integer :: N_MD_steps, N_loop
    integer :: n_extra_sorting
    integer :: n_threads
  
    type(threefry_rng_t), allocatable :: state(:)
  
    double precision :: skin
    double precision :: rsq
  
negative.code:   solvent% force_old = solvent% force
    colloids% force_old = colloids% force
  
    solvent% pos_old = solvent% pos
    colloids% pos_old = colloids% pos
    write(*,*) 'Running for', N_loop, 'loops'
    do i = 1, N_loop
       if (modulo(i, 100)==0) write(*, '(i09)', advance='no') i
       md_loop: do j = 1, N_MD_steps
          call md_pos(solvent, dt)
          do k = 1, colloids% Nmax
  