anchor.diff:    janus_io%species_info%mode = H5MD_FIXED
     call janus_io%init(hfile, 'janus', colloids)
   
  -  solvent_io%force_info%store = .false.
  -  solvent_io%id_info%store = .false.
  -  solvent_io%position_info%store = .true.
  -  solvent_io%position_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  -  solvent_io%position_info%step = N_loop*N_MD_steps
  -  solvent_io%position_info%time = N_loop*N_MD_steps*dt
  -  solvent_io%image_info%store = .true.
  -  solvent_io%image_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  -  solvent_io%image_info%step = N_loop*N_MD_steps
  -  solvent_io%image_info%time = N_loop*N_MD_steps*dt
  -  solvent_io%velocity_info%store = .true.
  -  solvent_io%velocity_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  -  solvent_io%velocity_info%step = N_loop*N_MD_steps
  -  solvent_io%velocity_info%time = N_loop*N_MD_steps*dt
  -  solvent_io%species_info%store = .true.
  -  solvent_io%species_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  -  solvent_io%species_info%step = N_loop*N_MD_steps
  -  solvent_io%species_info%time = N_loop*N_MD_steps*dt
  -  call solvent_io%init(hfile, 'solvent', solvent)
  +  if (do_solvent_io) then
  +     solvent_io%force_info%store = .false.
  +     solvent_io%id_info%store = .false.
  +     solvent_io%position_info%store = .true.
  +     solvent_io%position_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  +     solvent_io%position_info%step = N_loop*N_MD_steps
  +     solvent_io%position_info%time = N_loop*N_MD_steps*dt
  +     solvent_io%image_info%store = .true.
  +     solvent_io%image_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  +     solvent_io%image_info%step = N_loop*N_MD_steps
  +     solvent_io%image_info%time = N_loop*N_MD_steps*dt
  +     solvent_io%velocity_info%store = .true.
  +     solvent_io%velocity_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  +     solvent_io%velocity_info%step = N_loop*N_MD_steps
  +     solvent_io%velocity_info%time = N_loop*N_MD_steps*dt
  +     solvent_io%species_info%store = .true.
  +     solvent_io%species_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
  +     solvent_io%species_info%step = N_loop*N_MD_steps
  +     solvent_io%species_info%time = N_loop*N_MD_steps*dt
  +     call solvent_io%init(hfile, 'solvent', solvent)
  +  end if
   
     do k = 1, solvent%Nmax
        solvent% vel(1,k) = threefry_normal(state(1))*sqrt(T)
  
positive.code: 
    call h5gclose_f(fields_group, error)
    call janus_io%close()
    call solvent_io%close()
    call radial_hist_el%close()
    call n_solvent_el%close()
    if (do_quaternion) then
  
negative.code:   type(polar_fields_t) :: polar
    integer(HID_T) :: polar_id
  
    logical :: do_rattle, do_read_links, do_lennard_jones, do_elastic, do_quaternion, do_solvent_io
    integer, allocatable :: links(:,:)
    double precision, allocatable :: links_d(:)
    double precision :: link_treshold
  