It's a wikifact that if $f:=(1+\sqrt{5})/2$, then the vertices of a 600-cell centered at the origin of 4-space with edges of length $1/f$ can be given as follows:

The first step is produce a $120\times4$ matrix $U$ with these vectors its rows. The game is to do this without using a for-loop. We set up our background,

x = QQ['x'].0
K.<f> = NumberField(x^2-x-1)
RR4 = VectorSpace(K,4) 
H.<i,j,k> = QuaternionAlgebra(K,-1,-1)

So we work with quaternions over $\rats(f)$, where $f$ is a zero of $f^2-f-1$. Due to an unfortunate feature of Python, we must remember now never to use $i$ or $j$ or $k$ as an index in a list comprehension. If we evaluate [i^2 for i in [1..5]] the value of $i$ becomes $5$ once complete.

Now we create $U$.

CP3 = CartesianProduct([-1,1],[-1,1],[-1,1],[0])
CP4 = CartesianProduct([-1,1],[-1,1],[-1,1],[-1,1])
fv1 = RR4([1/2,f/2,(f-1)/2,0])
fv2 = RR4([(f-1)/2,1/2,f/2,0])
fv3 = RR4([f/2,(f-1)/2,1/2,0])
xx= [fv1.pairwise_product(RR4(sign_v))\
    for sign_v in CP3]\
    +[fv2.pairwise_product(RR4(sign_v))\
    for sign_v in CP3]\
    +[fv3.pairwise_product(RR4(sign_v))\
    for sign_v in CP3]
mm = Matrix(xx)
mma = mm.matrix_from_columns([1,0,3,2])
mmb = mm.matrix_from_columns([2,3,0,1])
mmc = mm.matrix_from_columns([3,2,1,0])
MM = block_matrix([mm,mma,mmb,mmc], nrows=4)
ww = RR4([1/2,1/2,1/2,1/2])
xxd =\ 
    [ww.pairwise_product(RR4(sign_v))\
    for sign_v in CP4]+[RR4([1,0,0,0]),\  
    RR4([-1,0,0,0]),RR4([0,1,0,0]),RR4([0,-1,0,0]),\ 
    RR4([0,0,1,0]), RR4([0,0,-1,0]),RR4([0,0,0,1]),\
    RR4([0,0,0,-1])]
U = MM.stack(Matrix(xxd))

The best we can say for the above code is that it works. Interestingly enough, the columns of $U$ are orthogonal since $U^TU =30I$.

U.transpose()*U

[30  0  0  0]
[ 0 30  0  0]
[ 0  0 30  0]
[ 0  0  0 30]

Two vectors are adjacent in the 1-skeleton if their inner product is $1/f$.

DC = Graph([[1..120], lambda i,j: U[i-1].inner_product(U[j-1]) == 1/2*f])
DC.is_regular(), DC.is_vertex_transitive(), DC.degree()[0], DC.diameter()

(True, True, 12, 5)

We compute the orbits of the vertex stabilizer of $\aut{DC}$.

DCgrp = DC.automorphism_group( partition = [[1],[2..120]])
map( len, DCgrp.orbits())

[12, 12, 20, 20, 12, 12, 1, 30, 1]

We convert the rows of $U$ to quaternions.

qu = [uu[0]+uu[1]*i+uu[2]*j +uu[3]*k for uu in U.rows()]

The reduced trace of a quaternion $q$ is $q+q^*$, in other words it is twice its real part. The multiplicative order of a quaternion is determined by its real trace. We compute the partition of $qu$ by reduced trace

def ls_dict( tuples):
    dc = {} 
    for pair in tuples:
        if dc.has_key( pair[0]):
            dc[ pair[0]].append(pair[1])
        else:
            dc[ pair[0]] = [pair[1]]
    return dc

tuples = [(qq.reduced_trace(),qq) for qq in qu]
dc = ls_dict( tuples)
[(kit, len(dc[kit])) for kit in dc.keys()]

[(0, 30), (1, 20), (2, 1), (-f, 12), (-f + 1, 12), 
(f - 1, 12), (f, 12), (-2, 1), (-1, 20)]

Comparing the number of quaternions with given reduced trace with the sizes of the orbits of the vertex stabilizer, suggests that the two partitions are equal. You should verify this. You should also verify that the multiplicative order of an element of qu is determined by its reduced trace. (Use qq.order().) In particular the quaternions with reduced trace $f$ or $1-f$ have order 10. The quaternions in qu form a multiplicative group isomorphic to $SL(2,5)$. You can access this group in Sage by

T = SL(2,5)
T

Special Linear Group of degree 2 over Finite Field of size 5

The 600-cell is a Cayley graph for this group with a connection set consisting of a conjugacy class of elements of order 10. There are two such conjugacy classes, one formed by the elements with reduced trace $f$ and the other consisting of the elements of reduced trace $1-f$. We show now that the Cayley graph with respect to the first of these conjugacy classes is isomorphic to $DC$, you should verify that the second is.

rtf = [it for it in qu if it.reduced_trace() == f]
conn = Set([ it^(-1)*rtf[0]*it for it in qu]).list()
CG = Graph( [qu, lambda q1, q2: q2*q1^(-1) in conn])

and now

CG.is_isomorphic( DC)

True

The quaternion \[ th = (-1/2) + 1/2*i + 1/2*j + 1/2*k \] has order three and its orbits on $DC$ are cocliques. The quotient over its orbits is a graph on 40 vertices which, like DC, is locally an icosahedron. We confirm the first claim:

nbhd = DC.subgraph( vertices=DC[1])
nbhd.is_isomorphic( graphs.IcosahedralGraph())

True