EDIT: Improved distance computation efficiency a lot by using DistanceMatrix:

With[{r = 33 Sqrt[3]}, With[ (* Find integer-valued points on the sphere. *) {pts = {x, y, z} /. Solve[Element[{x, y, z}, Sphere[{1, 3, 5}, r]], Integers]}, (* Compute distance matrix for the points. *) DistanceMatrix[pts, DistanceFunction -> SquaredEuclideanDistance] // (* Find indices of edges which may be a part of a regular tetrahedron in this circumsphere by length. *) Position[(8/3) r^2] // (* Construct graph edges with coordinates as vertices. *) Map[UndirectedEdge @@ pts[[#]] &]]] // (* Find all 4-cliques. *) FindClique[Graph[#], {4}, All] & 

$r = 33\sqrt{3}$ evaluates now in 24 milliseconds on average on my laptop.

It should be noted that FindClique returns maximal cliques. For exact 4-clique subgraphs one should normally use FindIsomorphicSubgraph but that is not necessary here since these graphs don't have larger than 4-cliques as one can have only four points at equal distances from each other in the three-dimensional Euclidean space.

---

EDIT: Converted the search into a graph clique problem:

With[{r = 33 Sqrt[3]}, Select[ Subsets[ (* Find integer-valued points on the sphere. *) {x, y, z} /. Solve[Element[{x, y, z}, Sphere[{1, 3, 5}, r]], Integers], {2}], (* Select edges which may be a part of a regular tetrahedron in this circumsphere by length. *) SquaredEuclideanDistance @@ # == (8/3) r^2 &] // (* Construct a graph and find all 4-cliques. *) MapApply@UndirectedEdge // FindClique[#, {4}, All] &] 

This is not the fastest solution but it's quite concise, and for $r = 33\sqrt{3}$ evaluates in 0.56 seconds on my laptop.

---

EDIT: A new, much more scalable approach here, the original solution under it.

With[{r = 33 Sqrt[3]}, Select[ Subsets[ (* Find integer-valued points on the sphere. *) {x, y, z} /. Solve[Element[{x, y, z}, Sphere[{1, 3, 5}, r]], Integers], {2}], (* Select edges which may be a part of a regular tetrahedron in this circumsphere by length. *) SquaredEuclideanDistance @@ # == Evaluate[(8/3) r^2] &]] // Select[ (* Select from pairs of such edges so that *) Flatten[#, 1] & /@ Subsets[#, {2}], (* resulting vertices form a regular tetrahedron. *) Equal @@ (SquaredEuclideanDistance @@@ Subsets[#, {2}]) &] & // (* Delete redundant reorderings. *) DeleteDuplicatesBy[Sort] 

... this computes 130 solutions for $r = 33\sqrt{3}$ in less than 10 seconds on my laptop while my old code below requires almost four minutes to do it.

Solution for $r = 99\sqrt{3}$ below.

---

Old answer:

Brute-force search with a small tweak (find valid equilateral triangles first) is probably the most efficient solution when the sphere radius is 15 (original question):

With[ {(* Integer-valued points on the sphere. *) onsphere = {x, y, z} /. Solve[Element[{x, y, z}, Sphere[{1, 3, 5}, 15]], Integers], (* Function which returns True if all points in the list are pairwise equidistant. *) eqdist = Equal @@ (SquaredEuclideanDistance @@@ Subsets[#, {2}]) &}, Table[ (* Return valid tetrahedra. *) If[eqdist[Append[i, j]], Append[i, j], Nothing], (* Equilateral triangles on the sphere. *) {i, Select[Subsets[onsphere, {3}], eqdist]}, (* All individual points to test along with the triangles. *) {j, onsphere}]] // (* Delete redundant reorderings. *) Flatten[#, 1] & // DeleteDuplicatesBy[Sort] (* {} *) 

Result is empty which means there are no such solutions.

With sphere radius of $5\sqrt{3}$ 14 solutions are found in 0.6 seconds, matching @cvgmt's result.

I used PowersRepresentations as Daniel Lichtblau did here 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 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}*) 

r = 5 Sqrt[3]

(*same old code as Version 1, just different r*) (*{0.028064, Null}*) 

---