Random coordinates generator

I want to generate random coordinates for spheres in a box geometry. I'm using while loop and i have 2 condition. First one is the distance of coordinates. General distance formula was used so that the coordinates do not overlap. Second one is the porosity. When porosity is less than 0.45 generating should stop. My code is working correctly but when i reduce porosity condition less than 0.80 the algorithm stucks. It cannot reach that porosity even after hours. How can I improve it to generate coordinates faster? Any suggestions are appreciated.

#dist = math.sqrt(((x2-x1)**2) + ((y2-y1)**2) + ((z2-z1)**2))
import math
import random
import numpy as np
import matplotlib.pyplot as plt

A = 0.04       # x border.
B = 0.04       # y border.
C = 0.125      # z border.
V_total = A*B*C # volume
r = 0.006       # min distance of spheres.
radius = 0.003  # radius of spheres.
wall_distance = 0.003

Porosity = 1.0
coordinates = np.array([])
while Porosity >= 0.90:
    # coordinates
    x = random.uniform(wall_distance, A-wall_distance)
    y = random.uniform(wall_distance, B-wall_distance)
    z = random.uniform(wall_distance, C-wall_distance)
    coord1 = (x,y,z)
    if coordinates.shape[0] == 0: # add first one without condition
        coordinates = np.array([coord1])
    else:
        coordinates = np.vstack((coordinates, coord1))
    # seperate x,y,z and convert list for control
    d_x = coordinates[:,0]
    x = d_x.tolist()
    d_y = coordinates[:,1]
    y = d_y.tolist()
    d_z = coordinates[:,2]
    z = d_z.tolist()
    for j in range(len(y)):
        for k in range(j+1, len(z)):
            dist = math.sqrt(((x[j]-x[k])**2) + ((y[j]-y[k])**2) + ((z[j]-z[k])**2))
            if dist <= r:
                coordinates = coordinates[:-1, :] # if distance is less than r, remove last coordinate
    # check porosity
    V_spheres = (4/3) * (np.pi) * (radius**3) * len(coordinates)
    V_void = V_total - V_spheres
    Porosity = V_void / V_total

print("Porosity: {}".format(Porosity))
print("Number of spheres: {}".format(len(coordinates)))
    
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_xlim([0, A])
ax.set_ylim([0, B])
ax.set_zlim([0, C])
ax.set_title('Coordinates for spheres')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
p = ax.scatter(coordinates[:,0], coordinates[:,1], coordinates[:,2])
fig.colorbar(p)
plt.show()