anchor.diff:      end do
     end subroutine refuel
   
  +  subroutine concentration_field
  +    double precision :: dimer_orient(3),x(3),y(3),z(3)
  +    double precision :: solvent_pos(3,solvent% Nmax)
  +    double precision :: dz,r,theta,x_pos,y_pos,z_pos
  +    integer :: o
  +    integer :: check
  +
  +    dz = solvent_cells% edges(3)/400.d0
  +    dimer_orient = colloids% pos(:,2) - colloids% pos(:,1)
  +    z = dimer_orient/sqrt(dot_product(dimer_orient,dimer_orient))
  +    x = (/0.d0, 1.d0, -dimer_orient(2)/dimer_orient(3)/)
  +    x = x/sqrt(dot_product(x,x))
  +    y = (/z(2)*x(3)-z(3)*x(2),z(3)*x(1)-z(1)*x(3),z(1)*x(2)-z(2)*x(1)/)
  +    conc_z = 0
  +
  +    do o = 1, solvent% Nmax
  +       solvent_pos(:,o) = solvent% pos(:,o) - colloids% pos(:,1)
  +       x_pos = dot_product(x,solvent_pos(:,o))
  +       y_pos = dot_product(y, solvent_pos(:,o))
  +       z_pos = dot_product(z, solvent_pos(:,o))
  +       solvent_pos(:,o) = (/x_pos,y_pos,z_pos/)
  +    end do
  +    do o = 1, solvent% Nmax
  +       r = sqrt(solvent_pos(1,o)**2 + solvent_pos(2,o)**2)
  +       theta = atan(solvent_pos(2,o)/solvent_pos(1,o))
  +       solvent_pos(1,o) = r
  +       solvent_pos(2,o) = theta
  +       solvent_pos(3,o) = solvent_pos(3,o)+colloids% pos(3,1)
  +       if (solvent% species(o)==2) then
  +          check = floor(solvent_pos(3,o)/dz)
  +          conc_z(check) = conc_z(check) + 1
  +       end if 
  +    end do
  +    colloid_pos(:,1) = 0
  +    colloid_pos(3,1) = colloids% pos(3,1)
  +    colloid_pos(:,2) = 0
  +    colloid_pos(3,2) = d + colloids% pos(3,1)
  +  end subroutine concentration_field
  +
   end program setup_single_dimer
  
positive.code:           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 do
  
    write(*,*) 'n extra sorting', n_extra_sorting
    
  
negative.code: 
    call PTkill(config)
    
    open(15,file ='dimerdata_chemConc.txt')
    open(16,file ='dimerdata_Conc.txt')
    
    colloids% species(1) = 1
    colloids% species(2) = 2
    colloids% vel = 0
  