anchor.diff:  
     call colloids% init(2,2, mass) !there will be 2 species of colloids
   
  -  call hfile%create(args%output_file, 'RMPCDMD::single_dimer_pbc', &
  -       'N/A', 'Pierre de Buyl')
     call thermo_data%init(hfile, n_buffer=50, step=N_MD_steps, time=N_MD_steps*dt)
   
  +  call h5gclose_f(params_group, error)
     call PTkill(config)
     
     colloids% species(1) = 1
  
positive.code:   integer :: n_threads
    type(h5md_file_t) :: hfile
    type(h5md_element_t) :: dummy_element
    integer(HID_T) :: fields_group
    integer(HID_T) :: box_group
    type(thermo_t) :: thermo_data
    double precision :: temperature, kin_e
  
negative.code:   call time_change%init('change')
  
    call h5open_f(error)
    call hfile%create(args%output_file, 'RMPCDMD::single_dimer_pbc', &
         'N/A', 'Pierre de Buyl')
    call h5gcreate_f(hfile%id, 'parameters', params_group, error)
  
    prob = PTread_d(config,'probability', loc=params_group)
    bulk_rmpcd = PTread_l(config, 'bulk_rmpcd', 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)
    d = PTread_d(config, 'd', loc=params_group)
    
    tau = PTread_d(config, 'tau', loc=params_group)
    N_MD_steps = PTread_i(config, 'N_MD', loc=params_group)
    dt = tau / N_MD_steps
    N_loop = PTread_i(config, 'N_loop', loc=params_group)
  
    sigma_C = PTread_d(config, 'sigma_C', loc=params_group)
    sigma_N = PTread_d(config, 'sigma_N', loc=params_group)
  
    ! solvent index first, colloid index second, in solvent_colloid_lj
    epsilon(:,1) = PTread_dvec(config, 'epsilon_C', 2, loc=params_group)
    epsilon(:,2) = PTread_dvec(config, 'epsilon_N', 2, loc=params_group)
  
    sigma(:,1) = sigma_C
    sigma(:,2) = sigma_N
  