anchor.diff:  
     call neigh% update_list(colloids, solvent, max_cut+skin, solvent_cells, solvent_colloid_lj)
   
  -  call h5gcreate_f(hfile%id, 'fields', fields_group, error)
  -  call h5gclose_f(fields_group, error)
  -
     e1 = compute_force(colloids, solvent, neigh, solvent_cells% edges, solvent_colloid_lj)
     e2 = compute_force_n2(colloids, solvent_cells% edges, colloid_lj)
     e3 = compute_bead_force()
  
positive.code:   type(h5md_element_t) :: n_solvent_el
    type(h5md_element_t) :: link_angle_el
  
    integer, parameter :: block_length = 8
    type(axial_correlator_t) :: axial_cf
    integer(HID_T) :: correlator_group
  
negative.code:      mass(3*(i-1)+1) = rho * sigma_N**3 * 4 * pi/3
       mass(3*(i-1)+2) = rho * sigma_E**3 * 4 * pi/3
       mass(3*(i-1)+3) = rho * sigma_N**3 * 4 * pi/3
    end do
  
  
    allocate(radial_hist(N_colloids))
    do i = 1, N_colloids
       call radial_hist(i)%init(min(sigma_E, sigma_N)/4, max(sigma_E, sigma_N)*3, 100, n_species=N_species)
    end do
  
    call solvent% init(N,2, system_name='solvent') ! there will be 2 species of solvent particles
  