Genetic Drift Simulation Using Python

Question

Hi I spent a few minutes modelling this genetic drift simulator (ie. low time-effort). It's a very basic programme that assumes 2 alleles in the population and assumes a simple random allocation of alleles to offspring. My code is below :

import random
import matplotlib.pyplot as plt


size = 10000
gen = 50
#p = 0.5

population_p = [random.choice([0,1,2]) for _ in range(size)]
population_q = [2 - value_p for value_p in population_p]
#population_1 = [random.choice([2, 1, 0]) for _ in range(size)]

print(sum(population_p) / (size*2))
#print(population_p)
#print(population_q)
#print(population_1)


p_freq_time = []
q_freq_time = []
sum_time = []

for generation in range(gen):
    p_freq = sum(population_p) / (size*2)
    q_freq = sum(population_q) / (size*2)
    sum_freq = p_freq + q_freq
    p_freq_time.append(p_freq)
    q_freq_time.append(q_freq)
    sum_time.append(sum_freq)
    
    
    next_generation = []
    for _ in range(size):
        parent_1 = random.choice(population_p)
        parent_2 = random.choice(population_p)
        low_bound = min(parent_1,parent_2)
        upper_bound = min(parent_1+parent_2,2)
        offspring = random.randint(low_bound,upper_bound)
        next_generation.append(offspring)   
    
    population_p = next_generation
    population_q = [2 - value_p for value_p in population_p]


plt.plot(range(gen), p_freq_time) #blue
#plt.plot(range(gen), q_freq_time) #orange
#plt.plot(range(gen), sum_time) #green
plt.xlabel('Generations')
plt.ylabel('p')
plt.ylim(-0.1,1.1)
plt.show()

Based on gene theory and the math on allele frequency, I am very sure that the theoretical framework is fine. However, is there a way to make the identification of offspring genetic code more random and reduce the chance of having a homozygous p slightly ? Because despite the randomisation, the fixation of p seems to happen very quickly even for large sample sizes.