Average Neighbouring Impurity Separation in a Random 1D chain 
I have a finite and discrete 1D chain (edit: linear chain, i.e. a straight line) of atoms, with unit separation, with a set number of impurities randomly distributed in the place of these atoms in the system. What I would like to do is describe the separation between neighbouring impurities (call it "D" which will always be an integer) statistically, and also to work out the average separation .
For example, the plot below was calculated from several thousand simulations of a chain of length 200 atoms and 10 impurities where the y-axis is the probability $P(D)$ of finding an impurity at distance D, and the x-axis is nearest impurity distance $D$. It kind of looks like a Poisson distributon, which one would expect since the system is discrete and random and a kind of counting exercise, but it doesn't work to well as a fit to the data points. It has been a long time since I did any statistics so I'm not sure how to start expressing what I found mathematically. Since I know the system length ($L = 200$) and the number of impurities ($N_i$) is a fair starting point the impurity density $\rho = N_i/L$ ? 
EDIT: The chain isn't allowed to self-intersect, it's a straight line in each case. The system I'm using above is a straight line of 200 evenly spaced atoms, and I'm distributing 10 impurities in the place of random atoms (e.g. at sites 4, 11, 54,...so there are still discrete steps between sites). The graph above is the result of finding the spacings between these impurity sites. 
EDIT 2: Attached a picture at the top

EDIT 3: Okay so it seems it could be my PRNG code causing problems. I'm using Fortran 95, here is the code:
`CALL RANDOM_SEED(size = n)
ALLOCATE(seed(n))

CALL SYSTEM_CLOCK(COUNT=clock) !!! intentionally slows it down to prevent succesive calls from returning the same number

do i = 1 , 1000000

end do
seed = clock + 37 * (/ (i - 1, i = 1, n) /)
CALL RANDOM_SEED(PUT = seed)

call random_number(x)`

EDIT 4!: Repeated this for a system size of 50 with 3 impurities (10,000 iterations), fitted a geometric distribution to it, as one can see immediately there is a huge variation. Is this explained by the fact for small n and L when we go from our first impurity (at site j) the probabilities to find the next impurity change drastically since effectively we are now looking for a system with n=2 and L= 50-j, which would be really sensitive to the location of the first impurity hence the huge variation. When we go to a huge system we can essentially treat the remainder of the chain as still being really long, and n and L are pretty much unchanged so this effect is masked and we see a geometric distribution arise? EDIT 5: Just did a bit of reading around Bernoulli Trials, is what I'm seeing and describing a binomial distribution hence the peak I see with smaller systems? Still that wouldn't explain why my 200 system is different from Ilmari's below...

 A: I'll conjecture that, at least in the limit as $L \to \infty$ while $\rho = N_i/L$ stays constant, $D$ will be geometrically distributed with parameter $\rho$.
This is because, for large $L$, we may essentially treat the states of each site as if they were i.i.d. Bernoulli random variables with probability $\rho$ of being an impurity.  Then the distribution of the distance from one impurity to the next (or, equivalently, from any arbitrary site to the next impurity) is the same as that for the waiting time before the first successful trial in a Bernoulli process with success rate $\rho$, which is what the geometric distribution describes.
To support this conjecture, here's my version of your plot:

The red crosses show the relative frequency of distances between nearest points over 1000 trials with $L=200$ and $N_i=10$, while the green line shows the expected frequencies given by $$P(D) = \rho (1-\rho)^{D-1}.$$
(Ps. From longer simulation runs, it seems evident that the effect of finite $L$ is to slightly increase the frequency of distances around the mean $1/\rho$ while reducing the frequency of longer distances above about $2/\rho$.  Intuitively, this seems perfectly reasonable, but even so, the deviation from the geometric distribution appears very slight, at least for these parameter values.)

For what it's worth, here's the Perl code I used to generate my (blipless) simulation data:

#!/usr/bin/perl
use warnings;
use strict;
use 5.010;

my $L = 200;         # number of sites
my $N = 10;          # number of impurities
my $reps = 1000;     # number of repeats

my %count;
for my $i (1 .. $reps) {
    # generate $N distinct random integers from 0 to $L-1 and sort them
    my %a;
    undef $a{int rand $L} while keys %a < $N;
    my @a = sort {$a <=> $b} keys %a;

    # count the differences between closest positions
    $count{ $a[$_] - $a[$_-1] }++ for 1 .. $#a;
}

say $_, "\t", $count{$_} / ($N-1) / $reps  for sort {$a <=> $b} keys %count;

A: I suspect you are seeing a flaw in linear congruential random number generators, that they generate successive random numbers restricted to various planes (see this infamous example). When you are generating 20 successive random integers, you are probably using the built in C rand() function, and the 20 points have correlations in their position.
Use the same program, but use a different generator, you can get one online. The problem should go away. You shouldn't use linear congruential generators for this sort of thing anyway.
