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.

I used PowersRepresentations as Daniel Lichtblau did here as it felt natural to me when I first saw this type of problem. With the new version below, I get the 130 solutions for $r=33~\sqrt{3}$ in less than 1 second on my laptop:

New version, r=33 Sqrt[3], a little simpler, a little faster, added some comments

New version, r = 5 Sqrt[3]

(*Same new code as in first codeblock, just different r*) (*{0.013104, Null}*) 

Old version, r = 33 Sqrt[3]

Old version, r = 5 Sqrt[3]

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

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}*) 

Version 2

An improvement on Version 1, added some comments

r= 33 Sqrt[3]

r = 5 Sqrt[3]

(*Same Version 2 code, just different r*) (*{0.013104, Null}*) 

Version 1

r = 33 Sqrt[3]

r = 5 Sqrt[3]

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

New version, r = 5 Sqrt[3]

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?

New version, r = 5 Sqrt[3]