I used `PowersRepresentations` as Daniel Lichtblau did [here](https://mathematica.stackexchange.com/a/289000/72953) as it felt natural to me when I first saw this type of problem.

With **Version 3**, I get the 130 solutions for $r=33~\sqrt{3}$ in less than 0.5 seconds on my laptop:


# Version 3 #


----------


Noticing that the `distances` matrix in **Version 2** was symmetrical, I roughly cut my time in half. I used Stelio's answer [here](https://mathematica.stackexchange.com/a/69236/72953) to create the `distances` matrix entries below the diagonal when I needed to index on it.

**r =33 Sqrt[3], using symmetry of distance matrix**

```
Clear["Global`*"];
AbsoluteTiming[center = {1, 3, 5};
 r = 33 Sqrt[3];
 (*get positive integer coordinates on r=15 sphere*)
 nnvals = PowersRepresentations[r^2, 3, 2];
 permVals = Flatten[Permutations /@ nnvals, 1];
 
 (*multiply the coords by all possible signs*)
 signs = Tuples[{-1, 1}, {3}];
 alltriples = Union[Flatten[Outer[Times, signs, permVals, 1], 1]];
 
 (*shift to center of sphere*)
 alltriples = # + center & /@ alltriples;
 
 (*side length from cvgmt's answer*)
 sideLength = r/(Sqrt[3/2]/2);
 
 (*note because of the +/- and permutations, lTrips is always even*)
 lTrips = Length@alltriples;
 (*take only the first half of alltriples, since the distance matrix \
is symmetrical*)
 upper = Take[alltriples, lTrips/2];
 (*calculate distances for first half of list to all of list*)
 distances = Outer[EuclideanDistance, upper, alltriples, 1];
 
 (*get list of coords that are sideLength away from each coord*)
 verticesWithoutSelf = (Flatten@Position[#, sideLength]) & /@ 
 distances;
 (*add the coord itself to its own list*)
 vertices = MapIndexed[Sort@Join[#2, #1] &, verticesWithoutSelf];
 
 (*get tetrahedron candidates,sort,and delete duplicates*)
 verts4 = Flatten[Sort@Subsets[#, {4}] & /@ vertices, 1];
 verts4 = DeleteDuplicates[verts4];
 
 (*create part of distance matrix below diagonal so we can index on \
it*)
 lowerDistances = 
 Reverse /@ (Transpose[Reverse /@ distances]) // Transpose;
 fullDistances = Join[distances, lowerDistances];
 
 (*grab indices from full distance matrix fullDistances*)
 distFun[v1_, v2_] := fullDistances[[v1, v2]];
 (*get distances between each vertex of the tetrahedron candidates*)
 edges = Outer[distFun[#1, #2] &, #, #, 1] & /@ verts4;
 (*get unique edge lengths*)
 edges = Sort[DeleteDuplicates[Flatten[#]]] & /@ edges;
 (*a regular tetrahedron will only have lengths sideLength and 0 \
(length to self)*)
 allSameSideLength = Flatten@Position[edges, {0, sideLength}];
 (*get tetrahedron vertex coordinates*)
 tetraVerts = verts4[[allSameSideLength]];
 tetrahedrons = Part[alltriples, #] & /@ tetraVerts;
]

(*{0.427283, Null}*)
```
**r = 5 Sqrt[3]**
```
(*same Version 3 code, just different r*)

(*{0.007682, Null}*)
```

As a small note, you don't actually have shift `alltriples` to the `center` of the sphere, you could apply this shift later to `tetrahedrons` and get the same thing since integer-integer = integer. This doesn't seem to make any difference in computation time however.

----------


# Version 2 #
An improvement on **Version 1**, added some comments

**r= 33 Sqrt[3]**
```
Clear["Global`*"];
AbsoluteTiming[
 center = {1, 3, 5};
 r = 33 Sqrt[3];
 (*get positive integer coordinates on r=15 sphere*)
 nnvals = PowersRepresentations[r^2, 3, 2];
 permVals = Flatten[Permutations /@ nnvals, 1];
 
 (*multiply the coords by all possible signs*)
 signs = Tuples[{-1, 1}, {3}];
 alltriples = Union[Flatten[Outer[Times, signs, permVals, 1], 1]];
 
 (*shift to center of sphere*)
 alltriples = # + center & /@ alltriples;
 
 (*side length from cvgmt's answer*)
 sideLength = r/(Sqrt[3/2]/2);
 
 (*calculate distance between all coords on the sphere*)
 distances = Outer[EuclideanDistance, alltriples, alltriples, 1];
 
 (*get list of coords that are sideLength away from each coord*)
 verticesWithoutSelf = (Flatten@Position[#, sideLength]) & /@ 
 distances;
 (*add the coord itself to its own list*)
 vertices = MapIndexed[Sort@Join[#2, #1] &, verticesWithoutSelf];
 
 (*get tetrahedron candidates, sort, and delete duplicates*)
 verts4 = Flatten[Sort@Subsets[#, {4}] & /@ vertices, 1];
 verts4 = DeleteDuplicates[verts4];
 
 (*I just made this function to get distances[[v1,v2]]*)
 distFun[v1_, v2_] := distances[[v1, v2]];
 (*get distances between each vertex of the tetrahedron candidates*)
 edges = Outer[distFun[#1, #2] &, #, #, 1] & /@ verts4;
 (*get unique edge lengths*)
 edges = Sort[DeleteDuplicates[Flatten[#]]] & /@ edges;
 
 (*a regular tetrahedron will only have lengths sideLength and 0 \
(length to self)*)
 allSameSideLength = Flatten@Position[edges, {0, sideLength}];
 (*get tetrahedron vertex coordinates*)
 tetraVerts = verts4[[allSameSideLength]];
 tetrahedrons = Part[alltriples, #] & /@ tetraVerts;]
Graphics3D[ConvexHullRegion /@ tetrahedrons]

(*{0.805232, Null}*)
``` 
For anyone wanting to further optimize for speed, the bottleneck is calculating the distances between all the integer coordinates:
```
AbsoluteTiming[
 distances = Outer[EuclideanDistance, alltriples, alltriples, 1];]
(*{0.762782, Null}*)
```
Also, the `AbsoluteTiming` of `distances` looks like it should be proportional to `Length[alltriples]^2`. Part of `alltriples` definition is a `Permutation`...this is probably not possible but, if we calculated the distances for one permutation of coordinates, could we draw conclusions about other permutations?

**r = 5 Sqrt[3]**
```
(*Same Version 2 code, just different r*)

(*{0.013104, Null}*)
```


----------
# Version 1 #

**r = 33 Sqrt[3]**

```
 Clear["Global`*"];
AbsoluteTiming[

 center = {1, 3, 5};
 r = 33 Sqrt[3];
 nnvals = PowersRepresentations[r^2, 3, 2];
 permVals = Flatten[Permutations /@ nnvals, 1];
 signs = Tuples[{-1, 1}, {3}];
 alltriples = Union[Flatten[Outer[Times, signs, permVals, 1], 1]];
 alltriples = # + center & /@ alltriples;
 
 sideLength = r/(Sqrt[3/2]/2);
 
 distances = Outer[EuclideanDistance, alltriples, alltriples, 1];
 
 tetras = Flatten[Table[
 candidateVertices = 
 Join[{i}, 
 Flatten[Position[distances[[i]], n_ /; n == sideLength]]];
 indices = Subsets[candidateVertices, {4}];
 Part[alltriples, #] & /@ indices
 , {i, distances // Length}], 1];
 sortedTetras = DeleteDuplicates[Sort /@ tetras];
 allSameSidePositions = 
 Flatten[Position[
 Table[Tally[
 DeleteCases[Flatten[Outer[EuclideanDistance, i, i, 1]], 
 0]][[All, 1]], {i, sortedTetras}], n_ /; n == {sideLength}]];
 tetrahedrons = sortedTetras[[allSameSidePositions]];
 ]
Graphics3D[ConvexHullRegion /@ tetrahedrons]
(*{2.0544, Null}*)
```
[![enter image description here][2]][2]


**r = 5 Sqrt[3]**

```
(*same old code as Version 1, just different r*)

(*{0.028064, Null}*)
```
[![enter image description here][1]][1]


----------


 [1]: https://i.sstatic.net/pGzn6.png
 [2]: https://i.sstatic.net/nqgIW.png