anchor.diff:  !! \param epsilon_N             interaction parameter of N sphere with both solvent species (2 elements)
   !! \param data_filename         filename for input Janus coordinates
   !! \param epsilon_colloid       interaction parameter for colloid-colloid interactions
  +!! \param reaction_radius       radius for the reaction around the Janus particle
   !! \param link_treshold         distance criterion for finding rigid-body links
   !! \param do_read_links         read link information from a file
   !! \param links_file            filename for the links data
  
positive.code:   solvent% force = 0
    solvent% species = 1
  
    call solvent_cells%init(L, 1.d0)
  
    call n_solvent_el%create_time(hfile%observables, 'n_solvent', &
         n_solvent, H5MD_LINEAR, step=N_MD_steps, &
         time=N_MD_steps*dt)
  
    call q_el%create_time(hfile%observables, 'q', &
         rigid_janus%q, H5MD_LINEAR, step=1, &
         time=dt)
  
    call janus_pos_el%create_time(hfile%observables, 'janus_pos', &
  
negative.code:            com_pos = sum((colloids%pos(:,:)+ &
                  colloids%image(:,:)*spread(solvent_cells%edges, dim=2, ncopies=colloids%Nmax) ) &
                  , dim=2)/ colloids%Nmax
             if (do_rattle .and. (.not. do_quaternion)) then
                rigid_janus%pos = com_pos
                rigid_janus%vel = v_com
             end if
             unit_r = get_unit_r()
             call axial_cf%add_fast((i-equilibration_loops)*N_MD_steps+j-1, v_com, unit_r)
             call janus_pos_el%append(rigid_janus%pos)
             call janus_vel_el%append(rigid_janus%vel)
  