anchor.diff:         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)
  +  if (do_quaternion) then
  +     call q_el%create_time(hfile%observables, 'q', &
  +          rigid_janus%q, H5MD_LINEAR, step=1, time=dt)
  +     call omega_body_el%create_time(hfile%observables, 'omega_body', &
  +          rigid_janus%omega_body, H5MD_LINEAR, step=1, time=dt)
  +     call omega_cf%init(block_length, get_n_blocks(block_length, 8, N_loop), dim=3)
  +  end if
   
     call u_el%create_time(hfile%observables, 'u', &
          unit_r, H5MD_LINEAR, step=1, &
  
positive.code:         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()
       end if
  
    end do
  
negative.code: 
    call h5gcreate_f(hfile%id, 'block_correlators', correlator_group, error)
    call axial_cf%write(correlator_group, N_MD_steps, N_MD_steps*dt, 1, dt)
    call write_correlator_block(correlator_group, 'omega_body_autocorrelation', omega_cf, N_MD_steps, N_MD_steps*dt)
    call h5gclose_f(correlator_group, error)
  
    ! write solvent coordinates for last step
  