Let us briefly talk about a probabilistic generalisation of k-means: the Gaussian Mixture Model (GMM).
In k-means, you carry out the following procedure:
- specify k centroids, initialising their coordinates randomly
- calculate the distance of each data point to each centroid
- assign each data point to its nearest centroid
- update the coordinates of the centroid to the mean of all points assigned to it
- iterate until convergence.
In a GMM, you carry out the following procedure:
- specify k multivariate Gaussians (termed components), initialising their mean and variance randomly
- calculate the probability of each data point being produced by each component (sometimes termed the responsibility each component takes for the data point)
- assign each data point to the component it belongs to with the highest probability
- update the mean and variance of the component to the mean and variance of all data points assigned to it
- iterate until convergence
You may notice the similarity between these two procedures. In fact, k-means is a GMM with fixed-variance components. Under a GMM, the probabilities (I think) you're looking for are the responsibilities each component takes for each data point.
There is a scikit-learn implementation of GMM available if you wanted to look into that, but I'm guessing you just want a quick way to amend your existing code, in which case, if you're happy to assume your clusters are fixed-variance Gaussians, you could transform your distance matrix element-wise as $y = e^{-x}$ (giving you an exponential fall-off), and then calculating the softmax over your columns (normalising your distribution so $P(Y=1) + P(Y=2) + ... + P(Y=k) = 1$).
It's worth pointing out that the assumption your clusters are fixed-variance Gaussians isn't necessarily valid. If your dimensions have wildly different scales, this may produce strange results, as dimensions with smaller-magnitude units will appear more "probable". Standardising your data before running your clustering procedure should remedy this.