anchor.diff:    call solvent_cells%init(L, 1.d0,has_walls = .true.)
     call vx% init(0.d0, solvent_cells% edges(3), L(3))
   
  -  allocate(rho_xy(N_species, L(2), L(1)))
  +  if (store_rho_xy) allocate(rho_xy(N_species, L(2), L(1)))
     call h5gcreate_f(hfile%id, 'fields', fields_group, error)
  -  call rho_xy_el%create_time(fields_group, 'rho_xy', rho_xy, ior(H5MD_LINEAR,H5MD_STORE_TIME), &
  +  if (store_rho_xy)  call rho_xy_el%create_time(fields_group, 'rho_xy', rho_xy, ior(H5MD_LINEAR,H5MD_STORE_TIME), &
          step=N_MD_steps, time=N_MD_steps*dt)
     call elem_vx% create_time(fields_group, 'vx', vx% data, ior(H5MD_TIME, H5MD_STORE_TIME))
     call elem_vx_count% create_time(fields_group, 'vx_count', vx% count, ior(H5MD_TIME, H5MD_STORE_TIME))
  
positive.code:   prob = PTread_d(config,'probability')
  
    number_of_angles = PTread_i(config, 'number_of_angles')
    dimer = PTread_l(config, 'dimer')
    if (dimer) then
       N_colloids = 2
  
negative.code: 
    call h5gclose_f(timers_group, error)
  
    if (store_rho_xy) call rho_xy_el%close()
    call elem_vx% close()
    call elem_vx_count% close()
    call dimer_io%close()
  