anchor.diff:    janus_io%id_info%store = .false.
     janus_io%position_info%store = .true.
     janus_io%position_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  -  janus_io%position_info%step = N_MD_steps
  -  janus_io%position_info%time = N_MD_steps*dt
  +  janus_io%position_info%step = colloid_sampling
  +  janus_io%position_info%time = colloid_sampling*dt
     janus_io%image_info%store = .true.
     janus_io%image_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  -  janus_io%image_info%step = N_MD_steps
  -  janus_io%image_info%time = N_MD_steps*dt
  +  janus_io%image_info%step = colloid_sampling
  +  janus_io%image_info%time = colloid_sampling*dt
     janus_io%velocity_info%store = .true.
     janus_io%velocity_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  -  janus_io%velocity_info%step = N_MD_steps
  -  janus_io%velocity_info%time = N_MD_steps*dt
  +  janus_io%velocity_info%step = colloid_sampling
  +  janus_io%velocity_info%time = colloid_sampling*dt
     janus_io%species_info%store = .true.
     janus_io%species_info%mode = H5MD_FIXED
     call janus_io%init(hfile, 'janus', colloids)
  +  call h5md_write_attribute(janus_io%position%id, 'head', head)
   
     solvent_io%force_info%store = .false.
     solvent_io%id_info%store = .false.
  
positive.code: program single_janus_pbc
    use rmpcdmd_module
    use hdf5
  
negative.code:   call timer_list%append(time_change)
  
    call h5open_f(error)
    call hfile%create(args%output_file, 'RMPCDMD::single_janus_pbc', &
         'N/A', 'Pierre de Buyl')
    call h5gcreate_f(hfile%id, 'parameters', params_group, error)
    call hdf5_util_write_dataset(params_group, 'seed', args%seed)
  
  
    prob = PTread_d(config,'probability', loc=params_group)
    bulk_rate = PTread_d(config,'bulk_rate', loc=params_group)
  
    L = PTread_ivec(config, 'L', 3, loc=params_group)
    rho = PTread_i(config, 'rho', loc=params_group)
    N = rho *L(1)*L(2)*L(3)
  
    T = PTread_d(config, 'T', 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)
    if (modulo(N_MD_steps, colloid_sampling) /= 0) then
       error stop 'colloid_sampling must divide N_MD with no remainder'
    end if
    dt = tau / N_MD_steps
    N_loop = PTread_i(config, 'N_loop', loc=params_group)
    equilibration_loops = PTread_i(config, 'equilibration_loops', loc=params_group)
  
    call thermo_data%init(hfile, n_buffer=50, step=N_MD_steps, time=N_MD_steps*dt)
  
    rattle_pos_tolerance = PTread_d(config, 'rattle_pos_tolerance', loc=params_group)
    rattle_vel_tolerance = PTread_d(config, 'rattle_vel_tolerance', loc=params_group)
  
    sigma = PTread_d(config, 'sigma_colloid', loc=params_group)
    sigma_v = sigma
    sigma_cut = sigma_v*2**(1.d0/6.d0)
    mass = rho * sigma**3 * 4 * 3.14159265/3
  
    epsilon = PTread_d(config, 'epsilon_colloid', loc=params_group)
  
    call colloid_lj% init(epsilon, sigma_v, sigma_cut)
  
    ! solvent index first, colloid index second, in solvent_colloid_lj
    sigma = PTread_d(config, 'sigma', loc=params_group)
    sigma_v = sigma
    sigma_cut = sigma_v*2**(1.d0/6.d0)
    max_cut = maxval(sigma_cut)
  
    epsilon(:,1) = PTread_dvec(config, 'epsilon_C', 2, loc=params_group)
    epsilon(:,2) = PTread_dvec(config, 'epsilon_N', 2, loc=params_group)
  
    call solvent_colloid_lj% init(epsilon, sigma_v, sigma_cut)
  
    call solvent% init(N, N_species, system_name='solvent') !there will be 2 species of solvent particles
  
    data_filename = PTread_s(config, 'data_filename', loc=params_group)
    link_treshold = PTread_d(config, 'link_treshold', loc=params_group)
  
    call h5gclose_f(params_group, error)
    call PTkill(config)
  
    call axial_cf%init(block_length, N_loop, N_loop*N_MD_steps)
  
    call colloids%init_from_file(data_filename, 'janus', H5MD_FIXED)
    head = get_head_idx(data_filename)
    colloids%image = 0
    colloids%vel = 0
    colloids%force = 0
  