I don't know if there is a polynomial-time algorithm (it feels like it might be NP-hard), but here are some plausible algorithmic approaches you might could consider, if you need to solve it in practice:
Heuristics
One well-explored algorithm is Furthest Point First (FPF). At each iteration, it chooses a point that is furthest from the set of points selected so far. Iterate $k$ times. As this is a greedy strategy, there is no reason to expect this to give an optimal answer or even close to optimal, and it was designed to optimize a slightly different objective function... but in some contexts it gives a reasonable approximation, so it could be worth a try.
FPF comes out of the literature on graph-based clustering and was introduced in the following research paper:
Teofilo F. Gonzalez. Clustering to minimize the maximum intercluster distance. Theoretical Computer Science, vol 38, pp.293-306, 1985.
You could try exploring the literature on graph-based clustering to see if anyone has studied your specific problem.
Exact algorithms
If you have this problem in practice and need an exact optimum solution, you could try to solve this using an ILP solver.
Here's how. Introduce 0-or-1 variables $x_i$, where $x_i$ indicates whether the $i$th vertex was selected, and 0-or-1 variables $y_{i,j}$, with the intended meaning that $y_{i,j}=1$ only if $x_i=1$ and $x_j=1$. Now maximize the objective function $\sum_{i,j} d(i,j) y_{i,j}$, subject to the constraints $\sum x_i \le k$ and $x_i \ge y_{i,j}$ and $x_j \ge y_{i,j}$. Now solve this ILP with an off-the-shelf ILP solver. As ILP is NP-hard, there's no guarantee this will be efficient, but it might work on some problem instances.
Another approach is to use weighted MAX-SAT. In particular, introduce boolean variables $x_i$, where $x_i$ is true if the $i$th vertex was selected, and variables $y_{i,j}$. The formula is $\phi \land \land_{i,j} y_{i,j}$, where $\phi$ must be true (its clauses have weight $W$ for some very large $W$) and each clause $y_{i,j}$ is given weight $d(i,j)$. Define the formula $\phi$ to be true if at most $k$ of the $x_i$'s are true (see here for details on how to do that) and if $y_{i,j}=x_i \land x_j$ for all $i,j$. Now the solution to this weighted MAX-SAT problem is the solution to the original problem, so you could try throwing a weighted MAX-SAT solver at the problem. The same caveats apply.
