anchor.diff:    write(*,*) ''
   
     write(*,*) 'n extra sorting', n_extra_sorting
  +
  +  ! create a group for block correlators and write the data
  +
  +  call h5gcreate_f(hfile%id, 'block_correlators', correlator_group, error)
  +
  +  call h5gcreate_f(correlator_group, 'mean_square_displacement', group, error)
  +  call h5md_write_dataset(group, 'value', msd%correlation)
  +  call h5md_write_dataset(group, 'count', msd%count)
  +  call h5md_write_dataset(group, 'step', N_MD_steps)
  +  call h5md_write_dataset(group, 'time', tau)
  +  call h5gclose_f(group, error)
  +
  +  call h5gcreate_f(correlator_group, 'orientation_autocorrelation', group, error)
  +  call h5md_write_dataset(group, 'value', oacf%correlation)
  +  call h5md_write_dataset(group, 'count', oacf%count)
  +  call h5md_write_dataset(group, 'step', N_MD_steps)
  +  call h5md_write_dataset(group, 'time', tau)
  +  call h5gclose_f(group, error)
  +
  +  call h5gclose_f(correlator_group, error)
  +
  +  ! write solvent coordinates for last step
   
     call solvent_io%position%append(solvent%pos)
     call solvent_io%velocity%append(solvent%vel)
  
positive.code:   type(thermo_t) :: thermo_data
    double precision :: temperature, kin_e
    double precision :: v_com(3)
    type(particle_system_io_t) :: dimer_io
    type(particle_system_io_t) :: solvent_io
    double precision :: bulk_rate
  
negative.code:      call dimer_io%velocity%append(colloids%vel)
       call dimer_io%image%append(colloids%image)
  
       com_pos = ( colloids%pos(:,1)+colloids%image(:,1)*solvent_cells%edges + &
            colloids%pos(:,2)+colloids%image(:,2)*solvent_cells%edges)
       call msd%add(i-1, correlate_block_distsq, xvec=com_pos)
  
       unit_r = rel_pos(colloids%pos(:,1), colloids%pos(:,2), solvent_cells%edges)
       unit_r = unit_r / norm2(unit_r)
       call oacf%add(i-1, correlate_block_dot, xvec=unit_r)
  
       solvent_cells% origin(1) = threefry_double(state(1)) - 1
       solvent_cells% origin(2) = threefry_double(state(1)) - 1
       solvent_cells% origin(3) = threefry_double(state(1)) - 1
  