Assessing Group Similarities and Dissimilarities Post PCA

To measure what features are the "drivers of difference between groups", you're going to need to frame this as a classification problem. Within this framework, you can use variable importance (and potentially coefficient values, depending on the model) to perform inference on drivers, i.e. to identify and rank them.

Applying PCA as a pre-processing step will make it significantly harder to identify what variables are drivers for specific classes. PCA is agnostic to class assignment: you can use PCA eigenvector components and loadings to interpret what variables are responsible for the bulk of the variance in your data, but that isn't actually what you want. Imagine if your classes were long elipsoidal clusters, each with the exact same covariance matrix but a different center (i.e. the elipses are "parellel"): the components given by PCA would be mainly influenced by the elipsoid shape (i.e. the shared covariance matrix) rather than by differences between groups. If you're interested in understanding the drivers of group differences, I strongly recommend you drop the PCA step.

I still don't have a good handle on exactly what you're hoping to get out of "which groups are most similar", but I suspect the "proximity" measure given by random forests would satisfy your need here. This measure is actually between observations, but you can take averages to get the expected proximity between groups. A benefit of using random forests here is that they have built-in variable importance measures, and you can even introspect them to understand the drivers behind individual observations.

Here's a little demo showing how to use a random forest model to detect drivers via variable importance, and measure group similarity via average inter-observation proximity:

First, set up the data and fit a model

library(randomForest)
data(iris)
set.seed(123)

unq_classes <- unique(iris$Species)
n_classes   <- length(unq_classes)

iris.rf <- randomForest(Species ~ ., data=iris, 
                        importance=TRUE,
                        proximity=TRUE)

Variable importances, rescaled to [0,1] with "1" indicating the most important variable:

var_imp <- importance(iris.rf)[,1:(n_classes+1)]
var_imp <- apply(var_imp, 2, function(m) m/max(m))

Here's the result (that last column is the marginalized importance):

                    setosa versicolor virginica MeanDecreaseAccuracy
Sepal.Length 0.2645863  0.2403256 0.2503047            0.3336813
Sepal.Width  0.1927240  0.0314708 0.1716495            0.1564093
Petal.Length 0.9525359  0.9589636 0.9356667            0.9549433
Petal.Width  1.0000000  1.0000000 1.0000000            1.0000000

And our mean proximities, again rescaled to [0,1] with 1 indicating the most similar pair of groups:

prx <- iris.rf$proximity
mean_prx <- matrix(NA, n_classes, n_classes)
for (i in 1:(n_classes-1)){
  for (j in (i+1):n_classes){
    cls_i <- iris$Species == unq_classes[i]
    cls_j <- iris$Species == unq_classes[j]
    mean_prx[j,i] <- mean(prx[cls_i, cls_j])
  }
}

mean_prx <- mean_prx/max(mean_prx, na.rm=TRUE)
rownames(mean_prx) <- unq_classes
colnames(mean_prx) <- unq_classes

Giving us:

                 setosa versicolor virginica
setosa               NA         NA        NA
versicolor 0.0267520374         NA        NA
virginica  0.0007778552          1        NA

Here's what the data looks like to put these results in context:

What does the distance between groups mean for your problem? Answering that problem would help you pick a distance metric.

Assuming you pick e.g. Euclidean distance, one simple way to find inter-group distance would be to first calculate the group centroid, i.e. average location of each group by averaging the locations of each group together. You can do this by just averaging all the 5x1 vectors belonging to a group, and repeat for each group. Then calculate the Euclidean distance between the resulting 6 centroids, giving you a 6x6 distance matrix.