anchor.diff:          v_com = sum(colloids%vel, dim=2)/colloids%Nmax
           call polar%update(com_pos, v_com, unit_r/norm2(unit_r), solvent, solvent_cells)
           call varia%tac()
  +
  +        if (do_quaternion) call omega_cf%add(i-equilibration_loops, correlate_block_dot, xvec=rigid_janus%omega_body)
  +
        end if
   
     end do
  
positive.code:        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 u_el%create_time(hfile%observables, 'u', &
         unit_r, H5MD_LINEAR, step=1, &
  
negative.code:   integer, dimension(N_species) :: n_solvent
    type(h5md_element_t) :: n_solvent_el
  
    type(h5md_element_t) :: q_el, omega_body_el, janus_pos_el, janus_vel_el, u_el
  
    type(histogram_t) :: radial_hist
    type(h5md_element_t) :: radial_hist_el
  