Oregon Curriculum Network
Discovering Math with Python

Quadrays and Graphene


By AlexanderAlUS - Own work, CC BY-SA 3.0, Link

"Graphene" refers to an hexagonal grid of cells, the vertexes being carbon atoms. However any hexagonal mesh, such as for game boards, might be referred to as a "graphene pattern".

Quadrays are explained in other Notebooks. Four basis vectors emanate to the corners of a volume 1 tetrahedron of edges 2R or 1D, in the canonical version, where R and D refer respectively to the Radius and Diameter of imaginary spheres packed together, giving this home base tetrahedron.

The Quadrays {2, 1, 1, 0}, meaning all 12 permutations of those numbers, fan out from (0,0,0,0) to the corners of a cuboctahedron.


In [1]:
from itertools import permutations
g = permutations((2,1,1,0))
unique = {p for p in g}  # set comprehension
print(unique)


{(0, 1, 1, 2), (1, 0, 1, 2), (2, 0, 1, 1), (0, 2, 1, 1), (0, 1, 2, 1), (1, 2, 1, 0), (1, 1, 2, 0), (2, 1, 1, 0), (1, 0, 2, 1), (1, 2, 0, 1), (2, 1, 0, 1), (1, 1, 0, 2)}

I have elsewhere used this fact to algorithmically generate consecutive shells of 12, 42, 92, 162... spheres (balls) respectively; a growing cuboctahedron of $10 S^{2} + 2$ balls per shell S = 1,2,3... (1 when S=0).

However suppose we don't want to grow the grid omni-directionally, but only in a plane. Each ball will be surrounded by six neighbors meaning at the center of a hexagon.

The cuboctahedron supplies four such hexagons i.e. its 24 edges comprise four hexagons orbiting the center. We may use any one of them.

The Algorithm

The algorithm begins with a planar subset of the vectors {2, 1, 1, 0} used to compute the six vertexes surrounding (0,0,0,0). We may call these six vertexes "carbons".

Then hop to neighboring hexagon centers (where no carbon is located) using an additional set of vectors.

From these new centers, compute the six surrounding carbons again, some of which will have already been found, as neighbors share fences, with three faces (centers) sharing each fence post (carbon).

Using the Python set object, the algorithm filters to keep only unique carbons.

Keep track of hexagon centers, a dual mesh, in a separate set.

(0,0,0,0) will be the first center (ring0).

If qrays r, s are 60 degrees apart on the same hexagon, pointing to neighboring carbons, then r + s will be the "hop" vector over the fence (edge) to the neighboring "yard" (face), or center.

Once we have six vertex vectors from a center, computing the six hop vectors (for jumping over the fences) will be a matter of summing pairs of adjacent (60 degree separated) vectors. We only keep new centers i.e. those of the next ring (see below).

What about edges?

As we go around a hexagon in 60 degree increments, say in a clockwise direction, we will be finding edges in terms of adjacent ball pairs.

To avoid redundancy, (ball_a, ball_b) -- any edge -- will be sorted. Any two quadrays may be ordered as 4-tuples e.g. (2, 1, 1, 0) is "greater than" (2, 1, 0, 1).

With unique representations of any edge, in the form of sorted tuples of qray namedtuples, we will be able to employ the same general technique employed with vertexes (carbons) and face centers: check the existing database for uniqueness and throw away (filter) anything already in the database. Sets will not allow duplicates.

The first step is to isolate six of the twelve from {2, 1, 1, 0} that define a hexagon.

One of the hexagons is TZOQXV. Do you see it in the above graphic? Another one is TYRQWS. If we regenerate all of the vectors A-Z mentioned above, we'll have a vocabulary suitable for graphene grid development, and then some.


In [2]:
from qrays import Qvector, IVM

A, B, C, D = Qvector((1,0,0,0)), Qvector((0,1,0,0)), Qvector((0,0,1,0)), Qvector((0,0,0,1))
E,F,G,H = B+C+D, A+C+D, A+B+D, A+B+C
I,J,K,L,M,N = A+B, A+C, A+D, B+C, B+D, C+D
O,P,Q,R,S,T = I+J, I+K, I+L, I+M, N+J, N+K; U,V,W,X,Y,Z = N+L, N+M, J+L, L+M, M+K, K+J

In [3]:
# two "beacons" of six spokes
hexrays = [T, Z, O, Q, X, V]              # to surrounding carbon atoms
hoprays = [T+Z, Z+O, O+Q, Q+X, X+V, V+T]  # to neighboring (vacant) hex centers

In [4]:
(T.angle(Z), Z.angle(O), O.angle(Q), Q.angle(X), X.angle(V), V.angle(T))


Out[4]:
(60.0, 60.0, 60.0, 60.0, 60.0, 60.0)

Lets verify that, going around the hexagon, each pair of consecutive hexrays is 60 degree apart. And ditto for hoprays, the vectors we'll use to jump over the fence to neighboring hexagon centers.


In [5]:
(hoprays[0].angle(hoprays[1]),
 hoprays[1].angle(hoprays[2]),
 hoprays[2].angle(hoprays[3]),
 hoprays[3].angle(hoprays[4]),
 hoprays[4].angle(hoprays[5]),
 hoprays[5].angle(hoprays[0]))


Out[5]:
(60.0, 60.0, 60.0, 60.0, 60.0, 60.0)

Looks like we're in business!

As with the growing cuboctahedron and the CCP packing, it makes sense to think in terms of consecutive rings.

The hexagonal coordination sequence is generated by:


In [6]:
def A008458(n):
    # OEIS number
    if n == 0:
        return 1
    return 6 * n
    
[A008458(x) for x in range(10)]


Out[6]:
[1, 6, 12, 18, 24, 30, 36, 42, 48, 54]

I will use this as a check as the algorithm generates multiple rings.


In [7]:
centers = {IVM(0,0,0,0)}  # center face
edges   = set()           # no duplicates permitted
carbons = set()

ring0 = [Qvector((0,0,0,0))]

def next_ring(ring):
    """
    Use only the most recently added hexagonal ring 
    of face centers to compute the next ring, moving 
    outward: 1, 6, 12, 18, 24...
    """
    new_faces = []
    for face in ring:
        verts = []
        
        # CARBONS
        for spoke in hexrays:
            v = face + spoke
            carbons.add(v.coords) # just the namedtuple is added to the set
            verts.append(v)
            
        # EDGES
        for bond in zip(verts, verts[1:] + [verts[0]]):
            # adding carbon-to-carbon bonds if not already in the set
            edge = tuple(sorted([bond[0].coords, bond[1].coords]))
            edges.add(edge)
            
        # CENTERS
        for jump in hoprays:
            neighbor = face + jump
            previous = len(centers)
            centers.add(neighbor.coords)
            if len(centers) > previous:  # if True, face is new
                new_faces.append(neighbor)
                
    return new_faces

In [8]:
def rings(n):
    prev = ring0
    for ring in range(n):
        print("Ring: {:3}  Number: {:4}".format(ring, len(prev)))
        nxt = next_ring(prev)
        prev = nxt

        
rings(12)


Ring:   0  Number:    1
Ring:   1  Number:    6
Ring:   2  Number:   12
Ring:   3  Number:   18
Ring:   4  Number:   24
Ring:   5  Number:   30
Ring:   6  Number:   36
Ring:   7  Number:   42
Ring:   8  Number:   48
Ring:   9  Number:   54
Ring:  10  Number:   60
Ring:  11  Number:   66

Note these are the expected numbers for consecutive rings.

Now that we have our database, it's time to generate some graphical output. As with the FCC, I'll use POV-Ray's scene description language and then render in POV-Ray. We just want to look at the edges and carbon atom vertexes.


In [9]:
sph = """sphere { %s 0.1 texture { pigment { color rgb <1,0,0> } } }"""
cyl = """cylinder { %s %s 0.05 texture { pigment { color rgb <1.0, 0.65, 0.0> } } }"""

def make_graphene(fname="../c6xty/graphene.pov", append=True):
    """
    Scan through carbons, edges, converting to XYZ and embedding
    in POV-Ray Scene Description Language
    """
    if append:
        pov = open(fname, "a")
    else:
        pov = open(fname, "w")

    # graphene will be included as a single object in the 
    # parent povray script, where lighting, camera position,
    # and background are defined
    
    print("#declare graphene = union{", file=pov)
    for atom in carbons:
        v = Qvector(atom).xyz()
        s = sph % "<{0.x}, {0.y}, {0.z}>".format(v)
        print(s, file=pov)
    for bond in edges:
        v0, v1 = bond
        v0 = Qvector(v0).xyz()
        v1 = Qvector(v1).xyz()
        c = cyl % ("<{0.x}, {0.y}, {0.z}>".format(v0), "<{0.x}, {0.y}, {0.z}>".format(v1))
        print(c, file=pov)
    print("}\n", file=pov)
        
make_graphene(append=False)

Success!

ADDENDUM

The graphrene grid needs twelve of its hexagons to become pentagons in order to lose the 720 degrees necessary to turn it into a concave / convex "global matrix" or "hexapent" grid.

This 720 degrees is characteristic of polyhedrons more generally and is sometimes referred to as "Descartes' Deficit". We may visualize the subtration of 720 degrees as the removal of one tetrahedron, thereby likewise begetting a tetrahedron (something with an inside and outside, not flat).