anchor.diff:       kin_e = (colloids% mass(1)*sum(colloids% vel(:,1)**2) + &
             colloids% mass(2)*sum(colloids% vel(:,2)**2))/2 + &
             sum(solvent% vel**2)/2
  -     write(15,*) colloids% pos + colloids% image * spread(solvent_cells% edges, dim=2, ncopies=colloids% Nmax), &
  -                 colloids% vel, e1+e2+kin_e
   
        call thermo_data%append(hfile, temperature, e1+e2, kin_e, e1+e2+kin_e)
   
        call dimer_io%position%append(colloids%pos)
  +     call dimer_io%velocity%append(colloids%vel)
  +     call dimer_io%image%append(colloids%image)
   
        solvent_cells% origin(1) = threefry_double(state(1)) - 1
        solvent_cells% origin(2) = threefry_double(state(1)) - 1
  
positive.code:   write(*,'(a16,f8.3)') main%name, main%total
  
  contains
  
    subroutine thermo_write
      write(*,'(1i16,6f16.3,1e16.8)') i, e1, e2, &
           kin_co, sum(solvent% vel**2)/2, &
           e1+e2+kin_co+sum(solvent% vel**2)/2, &
           compute_temperature(solvent, solvent_cells), &
           sqrt(dot_product(colloids% pos(:,1) - colloids% pos(:,2),colloids% pos(:,1) - colloids% pos(:,2))) - d
    end subroutine thermo_write
  
  
    subroutine flag_particles_nl
      double precision :: dist_to_C_sq
  
negative.code:   colloids% pos(:,1) = solvent_cells% edges/2.0
    colloids% pos(:,2) = solvent_cells% edges/2.0 
    colloids% pos(1,2) = colloids% pos(1,2) + d
  
    call h5gcreate_f(dimer_io%group, 'box', box_group, error)
    call h5md_write_attribute(box_group, 'dimension', 3)
    call dummy_element%create_fixed(box_group, 'edges', solvent_cells%edges)
    call h5gclose_f(box_group, error)
  
    call solvent% random_placement(solvent_cells% edges, colloids, solvent_colloid_lj)
  
    call solvent% sort(solvent_cells)
  