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 dummy_element%create_fixed(box_group, 'edges', solvent_cells%edges)
    call h5gclose_f(box_group, error)
  
    call h5gcreate_f(solvent_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)
  
    call h5gcreate_f(hfile%id, 'fields', fields_group, error)
    call radial_hist%init(0.d0, 5*sigma, 100, 2)
  
negative.code:   ! write solvent coordinates for last step
  
    call thermo_data%append(hfile, temperature, e1+e2+e3, kin_e, e1+e2+e3+kin_e, v_com, add=.false., force=.true.)
    if (do_solvent_io) then
       call solvent_io%position%append(solvent%pos)
       call solvent_io%velocity%append(solvent%vel)
       call solvent_io%image%append(solvent%image)
       call solvent_io%species%append(solvent%species)
       call solvent_io%close()
    end if
  
    ! store polar fields
  
  