anchor.diff:         z_pos = dot_product(z, solvent_pos(:,o))
          solvent_pos(:,o) = (/x_pos,y_pos,z_pos/)
       end do
  +    !$omp parallel do private(r, theta, check)
       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))
  
positive.code:   double precision :: x(3)
    
    call random_number(rndnumbers)
    
    do  r = 1,solvent% Nmax
       if (solvent% species(r) == 1) then
         x = rel_pos(colloids% pos(:,1),solvent% pos(:,r),solvent_cells% edges) 
  
negative.code:        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)
            !$omp atomic
            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
  