anchor.diff:       T = T_init + (i-1)*(T_final-T_init)/(N_loop-1)
        call mpcd_at_step(solvent, solvent_cells, state, T)
   
  -     if (modulo(i,100)==0) &
  -     write(15,*) colloids% pos + colloids% image * spread(solvent_cells% edges, dim=2, ncopies=colloids% Nmax)
  -     if (modulo(i,10)==0) &
  -     write(*,'(6f16.3)') e1, e2, mass*sum(colloids% vel**2)/2, sum(solvent% vel**2)/2, &
  -         e1+e2+mass*sum(colloids% vel**2)/2+sum(solvent% vel**2)/2, &
  -         compute_temperature(solvent, solvent_cells)
  +     if (modulo(i,colloid_sampling)==0) then
  +        call colloids_io%position%append(colloids%pos)
  +        call colloids_io%image%append(colloids%image)
  +     end if
  +
  +     temperature = compute_temperature(solvent, solvent_cells)
  +     total_mass = 0
  +     kin_e = sum(solvent% vel**2)/2
  +     v_com = sum(solvent% vel, dim=2)
  +     do k = 1, colloids% Nmax
  +        j = colloids%species(k)
  +        if (j==0) cycle
  +        local_mass = colloids%mass(j)
  +        kin_e = kin_e + local_mass*sum(colloids% vel(:,k)**2)/2
  +        v_com = v_com + local_mass*colloids%vel(:,k)
  +        total_mass = total_mass + local_mass
  +     end do
  +     v_com = v_com / (solvent%Nmax + total_mass)
  +
  +     call thermo_data%append(hfile, temperature, e1+e2, kin_e, e1+e2+kin_e, v_com)
   
     end do
   
  +  call thermo_data%append(hfile, temperature, e1+e2, kin_e, e1+e2+kin_e, &
  +       v_com, add=.false., force=.true.)
  +
     write(*,*) 'n extra sorting', n_extra_sorting
   
  +  call colloids_io%close()
  +  call hfile%close()
     call h5close_f(error)
   
   end program n_colloids_pbc
  
positive.code:   double precision :: so_max, co_max
  
    double precision :: e1, e2
    double precision :: tau, dt, T, T_init, T_final
    integer :: N_MD_steps, N_loop
    integer :: n_extra_sorting
    integer :: n_threads
  
    type(threefry_rng_t), allocatable :: state(:)
  
    double precision :: skin
    double precision :: rsq
  
negative.code:   call threefry_rng_init(state, args%seed)
  
    call h5open_f(error)
    call hfile%create(args%output_file, 'RMPCDMD::n_colloids_pbc', &
         RMPCDMD_REVISION, 'Pierre de Buyl')
    call h5gcreate_f(hfile%id, 'parameters', params_group, error)
    call hdf5_util_write_dataset(params_group, 'seed', args%seed)
  
    L = PTread_ivec(config, 'L', 3, loc=params_group)
    rho = PTread_i(config, 'rho', loc=params_group)
    T_init = PTread_d(config, 'T', loc=params_group)
    T = T_init
    T_final = PTread_d(config, 'T_final', loc=params_group)
    tau = PTread_d(config, 'tau', loc=params_group)
    N_MD_steps = PTread_i(config, 'N_MD', loc=params_group)
    colloid_sampling = PTread_i(config, 'colloid_sampling', loc=params_group)
    dt = tau / N_MD_steps
    N_loop = PTread_i(config, 'N_loop', loc=params_group)
  
    N_colloids = PTread_i(config, 'N_colloids', loc=params_group)
    epsilon = PTread_d(config, 'epsilon', loc=params_group)
    sigma = PTread_d(config, 'sigma', loc=params_group)
    sigma_cut = sigma*2.d0**(1.d0/6.d0)
  
    call thermo_data%init(hfile, n_buffer=50, step=N_MD_steps, time=N_MD_steps*dt)
  
    mass = rho * sigma**3 * 4 * 3.141/3
    call colloids%init(N_colloids, 1, [mass])
    colloids%species = 1
  
    colloids_io%force_info%store = .false.
    colloids_io%id_info%store = .false.
    colloids_io%velocity_info%store = .false.
    colloids_io%position_info%store = .true.
    colloids_io%position_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
    colloids_io%position_info%step = N_MD_steps
    colloids_io%position_info%time = N_MD_steps*dt
    colloids_io%image_info%store = .true.
    colloids_io%image_info%mode = ior(H5MD_LINEAR,H5MD_STORE_TIME)
    colloids_io%image_info%step = N_MD_steps
    colloids_io%image_info%time = N_MD_steps*dt
    colloids_io%species_info%store = .true.
    colloids_io%species_info%mode = H5MD_FIXED
    call colloids_io%init(hfile, 'colloids', colloids)
  
    call h5gcreate_f(colloids_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)
  
    N = rho*L(1)*L(2)*L(3) - int(rho*4*3.142/3 * sigma**3*colloids%Nmax)
  
    write(*,*) N, 'solvent particles'
  
    if (N <= 0) error stop 'Not enough volume available for solvent'
  
    call solvent_colloid_lj% init( reshape( [ epsilon ], [1, 1] ), &
         reshape( [ sigma ], [1, 1] ), reshape( [ sigma_cut ], [1, 1] ) )
  
    epsilon = PTread_d(config, 'epsilon_colloids', loc=params_group)
    sigma_cut = 3*sigma
  
    call h5gclose_f(params_group, error)
    call PTkill(config)
  
    call colloid_lj% init( reshape( [ epsilon ], [1, 1] ), &
         reshape( [ 2*sigma ], [1, 1] ), reshape( [ 2*sigma_cut ], [1, 1] ) )
  