Simulating substitution rate of neutral mutations
I don't fully understand how your model works

I don't fully understand how you model your population. Is it a standard $oop$ (object-oriented programming) where you simulate each individual? Or is it a simulation where you already use some mathematical model? You could eventually copy-paste your code. it is probably not very long, right? What language did you use? Do you have only one locus? Do you simulate several positions within these locus? Do you allow for back mutations?

It seems that you allow the population size to change over time and I don't quite understand how your parameter $K$ does limit the growth rate though. You should write down some equations about that. In any case, change in population size will modify the probability of mutation fixation.

Here is probably some important modifications to do on your code

It is a bit hard to comment on your simulations without actually seeing some equations or seeing the code but I really think that you'll need two important modifications in your code.


fix population size $N$. 
To simulate reproduction, randomly choose (uniform distribution) individuals (with replacement) from the population. The chosen individuals form the new generation. For example:


In R:

new_pop = sample(old_pop, N, replace=T)


or directly

pop = sample(pop, N, replace=T)


In Python:

from random import randrange
new_pop = []
for i in xrange(N):
    new_pop.append(old_pop[randrange(N)])


When you get your results, don't forget that you're talking about a stochastic process and you may want to run enough simulation or enough independent loci in order to be able to compare with the predictions from mathematical models.