Is the function $k(g,h) = \frac{1}{1+\lvert gh^{-1}\rvert}$ positive definite? Let $G$ be a finite group, $S \subset G$ a generating set, closed under taking inverses, and $\lvert\cdot\rvert$ the word length with respect to this set $S$.

Question. Is the function $k(g,h) = \frac{1}{1+\lvert gh^{-1}\rvert}$ positive definite, for $g,h \in G$?

A positive answer would allow every finite group $G$ to be embedded as a metric space in Euclidean space with $d(g,h) = \sqrt{2-2k(g,h)}$, and this could be useful in Geometric Group Theory.
Motivation:
This matrix $k(g,h)$ plays a role in a group theoretic reformulation of the Lagarias inequality:
What properties characterize the function $L(x) = x+\exp(x) \log(x)$?
Also asked on MSE, since I was not sure if it is research relevant:
https://math.stackexchange.com/questions/4126491/is-the-function-kg-h-frac11gh-1-positive-definite
 A: [NB. Addendum below gives a counterexample to OP's question.]
Computationally, this appears to be true for $S_n$, $2 \le n \le 6$.
The MATLAB code I used to check this is below (using the case $n = 4$ to eyeball the Cayley graph). By Cayley's theorem and appropriately tweaking the code (essentially, the arrays gen and G), you should be able to check other small groups with this as well (this is just exploiting the permutation representation, albeit very expensively).
As an aside, it also appears that the kernel $\exp(-\beta|gh^{-1}|)$ is positive definite on $S_n$ for all $\beta > 0$ using adjacent transpositions, but the story is more delicate for all transpositions (see links in reference).
n = 6;

% Generating set of S_n: adjacent tranpositions
gen = ones(n-1,1)*(1:n);
gen(1:n:(n-1)^2) = 2:n;
gen(n:n:((n-1)*n)) = 1:(n-1);

% S_n itself
G = perms(1:n);

% Indices of generating set
genInd = find(ismember(G,gen,'rows'));

% Represent symmetric group using permutation matrices
P = cell(size(G,1),1);
for j = 1:numel(P)
    P{j} = zeros(n);
    for k = 1:n
        P{j}(k,G(j,k)) = 1;
    end
end

% Multiplication table in terms of group indices
mul = nan(size(G,1));
for j = 1:numel(P)
    for k = 1:numel(P)
        [prod_jk,~] = find(P{j}*P{k});
        mul(j,k) = find(ismember(G,prod_jk','rows'));
    end
end

% Get indices of inverses from this
[ver,~] = find(mul==find(ismember(G,1:n,'rows')));  % don't overload inv

% Division table
div = mul(:,ver);

% Adjacency matrix of Cayley graph
A = ismember(div,genInd);

% Cayley graph (Bruhat for this particular case)
nodenames = join([repmat("x",[size(G,1),1]),string(G)],'');
cayleyGraph = digraph(A,nodenames);

% Distances equal word lengths
d = distances(cayleyGraph);

% Form kernel
ker = 1./(1+d);

% Inspect spectrum visually to check positive definiteness
spec = eig(ker)'

% Optional: plot Cayley graph for assurance
% figure;
% plot(cayleyGraph);

ADDENDUM:
For $S_5$, using all transpositions shows this is false. The key bit to change in the above code is to put this generating set in after G is defined:
% All transpositions
transInd = find(sum(G~=ones(size(G,1),1)*(1:n),2)==2);
gen = G(transInd,:);

You can verify that gen obeys the hypotheses by checking that transInd and ver(transInd) have the same indices. In this case, spec shows an eigenvalue of -0.0333.
