GrapheneWithQrays.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Oregon Curriculum Network](http://www.4dsolutions.net/ocn) <br />\n",
    "[Discovering Math with Python](Introduction.ipynb)\n",
    "\n",
    "# Quadrays and Graphene\n",
    "\n",
    "<p><a href=\"https://commons.wikimedia.org/wiki/File:Graphen.jpg#/media/File:Graphen.jpg\"><img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/9/9e/Graphen.jpg/1200px-Graphen.jpg\" alt=\"Graphen.jpg\"></a><br>By <a href=\"//commons.wikimedia.org/w/index.php?title=User:AlexanderAlUS&amp;action=edit&amp;redlink=1\" class=\"new\" title=\"User:AlexanderAlUS (page does not exist)\">AlexanderAlUS</a> - <span class=\"int-own-work\" lang=\"en\">Own work</span>, <a href=\"https://creativecommons.org/licenses/by-sa/3.0\" title=\"Creative Commons Attribution-Share Alike 3.0\">CC BY-SA 3.0</a>, <a href=\"https://commons.wikimedia.org/w/index.php?curid=11294534\">Link</a></p>\n",
    "\n",
    "\"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\".\n",
    "\n",
    "Quadrays are explained [in other Notebooks](QuadraysJN.ipynb).  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.\n",
    "\n",
    "![quadrays](https://upload.wikimedia.org/wikipedia/commons/9/99/Quadray.gif)\n",
    "\n",
    "<a data-flickr-embed=\"true\"  href=\"https://farm5.staticflickr.com/4143/4949799682_327b33e8d5.jpg\" title=\"Tetrahedron\"><img src=\"https://farm5.staticflickr.com/4143/4949799682_327b33e8d5.jpg\" width=\"375\" height=\"500\" alt=\"Tetrahedron\"></a><script async src=\"//embedr.flickr.com/assets/client-code.js\" charset=\"utf-8\"></script>\n",
    "\n",
    "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.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{(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)}\n"
     ]
    }
   ],
   "source": [
    "from itertools import permutations\n",
    "g = permutations((2,1,1,0))\n",
    "unique = {p for p in g}  # set comprehension\n",
    "print(unique)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I have [elsewhere](Generating%20the%20FCC.ipynb) 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).\n",
    "\n",
    "![Image of Cubocta](http://www.4dsolutions.net/ocn/graphics/cubanim.gif)\n",
    "\n",
    "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.\n",
    "\n",
    "The cuboctahedron supplies four such hexagons i.e. its 24 edges comprise four hexagons orbiting the center.  We may use any one of them.\n",
    "\n",
    "## The Algorithm\n",
    "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\".\n",
    "\n",
    "Then hop to neighboring hexagon centers (where no carbon is located) using an additional set of vectors.  \n",
    "\n",
    "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).  \n",
    "\n",
    "Using the Python set object, the algorithm filters to keep only unique carbons.  \n",
    "\n",
    "Keep track of hexagon centers, a dual mesh, in a separate set.  \n",
    "\n",
    "(0,0,0,0) will be the first center (ring0).\n",
    "\n",
    "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.  \n",
    "\n",
    "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).\n",
    "\n",
    "What about edges?\n",
    "\n",
    "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.  \n",
    "\n",
    "To avoid redundancy, (ball_a, ball_b) -- any edge -- [will be sorted](https://github.com/4dsolutions/SAISOFT/blob/master/OrderingPolys.ipynb).  Any two quadrays may be ordered as 4-tuples e.g. (2, 1, 1, 0) is \"greater than\" (2, 1, 0, 1).  \n",
    "\n",
    "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.\n",
    "\n",
    "The first step is to isolate six of the twelve from {2, 1, 1, 0} that define a hexagon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<table class=\"multicol\" role=\"presentation\" style=\"border-collapse: collapse; padding: 0; border: 0; background:transparent; width:100%;\"><tbody><tr>\n",
    "<td style=\"text-align: left; vertical-align: top;\">\n",
    "<table class=\"wikitable\">\n",
    "\n",
    "<tbody><tr>\n",
    "<th>Shape\n",
    "</th>\n",
    "<th>Volume\n",
    "</th>\n",
    "<th>Vertex Inventory (sum of Quadrays)\n",
    "</th></tr>\n",
    "<tr>\n",
    "<td>Tetrahedron\n",
    "</td>\n",
    "<td>1\n",
    "</td>\n",
    "<td>A,B,C,D\n",
    "</td></tr>\n",
    "<tr>\n",
    "<td>Inverse Tetrahedron\n",
    "</td>\n",
    "<td>1\n",
    "</td>\n",
    "<td>E,F,G,H = B+C+D, A+C+D, A+B+D, A+B+C\n",
    "</td></tr>\n",
    "<tr>\n",
    "<td>Duo-Tet Cube\n",
    "</td>\n",
    "<td>3\n",
    "</td>\n",
    "<td>A through H\n",
    "</td></tr>\n",
    "<tr>\n",
    "<td>Octahedron\n",
    "</td>\n",
    "<td>4\n",
    "</td>\n",
    "<td>I,J,K,L,M,N = A+B, A+C, A+D, B+C, B+D, C+D\n",
    "</td></tr>\n",
    "<tr>\n",
    "<td>Rhombic Dodecahedron\n",
    "</td>\n",
    "<td>6\n",
    "</td>\n",
    "<td>A through N\n",
    "</td></tr>\n",
    "<tr>\n",
    "<td>Cuboctahedron\n",
    "</td>\n",
    "<td>20\n",
    "</td>\n",
    "<td>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\n",
    "</td></tr></tbody></table>\n",
    "<p>&#32;\n",
    "</p>\n",
    "</td>\n",
    "<td style=\"text-align: left; vertical-align: top;\">\n",
    "<p><a href=\"https://upload.wikimedia.org/wikipedia/commons/d/dc/Povlabels.gif\" class=\"image\" title=\"Points A-Z\"><img alt=\"Points A-Z\" src=\"https://upload.wikimedia.org/wikipedia/commons/d/dc/Povlabels.gif\" width=\"320\" height=\"240\" data-file-width=\"320\" data-file-height=\"240\" /></a>\n",
    "&#32;\n",
    "</p>\n",
    "</td></tr></tbody></table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from qrays import Qvector, IVM\n",
    "\n",
    "A, B, C, D = Qvector((1,0,0,0)), Qvector((0,1,0,0)), Qvector((0,0,1,0)), Qvector((0,0,0,1))\n",
    "E,F,G,H = B+C+D, A+C+D, A+B+D, A+B+C\n",
    "I,J,K,L,M,N = A+B, A+C, A+D, B+C, B+D, C+D\n",
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# two \"beacons\" of six spokes\n",
    "hexrays = [T, Z, O, Q, X, V]              # to surrounding carbon atoms\n",
    "hoprays = [T+Z, Z+O, O+Q, Q+X, X+V, V+T]  # to neighboring (vacant) hex centers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(60.0, 60.0, 60.0, 60.0, 60.0, 60.0)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(T.angle(Z), Z.angle(O), O.angle(Q), Q.angle(X), X.angle(V), V.angle(T))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(60.0, 60.0, 60.0, 60.0, 60.0, 60.0)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(hoprays[0].angle(hoprays[1]),\n",
    " hoprays[1].angle(hoprays[2]),\n",
    " hoprays[2].angle(hoprays[3]),\n",
    " hoprays[3].angle(hoprays[4]),\n",
    " hoprays[4].angle(hoprays[5]),\n",
    " hoprays[5].angle(hoprays[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks like we're in business!\n",
    "\n",
    "As with the growing cuboctahedron and the CCP packing, it makes sense to think in terms of consecutive rings.\n",
    "\n",
    "The [hexagonal coordination sequence](https://oeis.org/A008458) is generated by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 6, 12, 18, 24, 30, 36, 42, 48, 54]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def A008458(n):\n",
    "    # OEIS number\n",
    "    if n == 0:\n",
    "        return 1\n",
    "    return 6 * n\n",
    "    \n",
    "[A008458(x) for x in range(10)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I will use this as a check as the algorithm generates multiple rings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "centers = {IVM(0,0,0,0)}  # center face\n",
    "edges   = set()           # no duplicates permitted\n",
    "carbons = set()\n",
    "\n",
    "ring0 = [Qvector((0,0,0,0))]\n",
    "\n",
    "def next_ring(ring):\n",
    "    \"\"\"\n",
    "    Use only the most recently added hexagonal ring \n",
    "    of face centers to compute the next ring, moving \n",
    "    outward: 1, 6, 12, 18, 24...\n",
    "    \"\"\"\n",
    "    new_faces = []\n",
    "    for face in ring:\n",
    "        verts = []\n",
    "        \n",
    "        # CARBONS\n",
    "        for spoke in hexrays:\n",
    "            v = face + spoke\n",
    "            carbons.add(v.coords) # just the namedtuple is added to the set\n",
    "            verts.append(v)\n",
    "            \n",
    "        # EDGES\n",
    "        for bond in zip(verts, verts[1:] + [verts[0]]):\n",
    "            # adding carbon-to-carbon bonds if not already in the set\n",
    "            edge = tuple(sorted([bond[0].coords, bond[1].coords]))\n",
    "            edges.add(edge)\n",
    "            \n",
    "        # CENTERS\n",
    "        for jump in hoprays:\n",
    "            neighbor = face + jump\n",
    "            previous = len(centers)\n",
    "            centers.add(neighbor.coords)\n",
    "            if len(centers) > previous:  # if True, face is new\n",
    "                new_faces.append(neighbor)\n",
    "                \n",
    "    return new_faces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ring:   0  Number:    1\n",
      "Ring:   1  Number:    6\n",
      "Ring:   2  Number:   12\n",
      "Ring:   3  Number:   18\n",
      "Ring:   4  Number:   24\n",
      "Ring:   5  Number:   30\n",
      "Ring:   6  Number:   36\n",
      "Ring:   7  Number:   42\n",
      "Ring:   8  Number:   48\n",
      "Ring:   9  Number:   54\n",
      "Ring:  10  Number:   60\n",
      "Ring:  11  Number:   66\n"
     ]
    }
   ],
   "source": [
    "def rings(n):\n",
    "    prev = ring0\n",
    "    for ring in range(n):\n",
    "        print(\"Ring: {:3}  Number: {:4}\".format(ring, len(prev)))\n",
    "        nxt = next_ring(prev)\n",
    "        prev = nxt\n",
    "\n",
    "        \n",
    "rings(12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a data-flickr-embed=\"true\"  href=\"https://www.flickr.com/photos/kirbyurner/2949087863/in/album-72157612943105800/\" title=\"Global Matrix\"><img src=\"https://farm4.staticflickr.com/3249/2949087863_aea4cecb7f.jpg\" width=\"500\" height=\"375\" alt=\"Global Matrix\"></a><script async src=\"//embedr.flickr.com/assets/client-code.js\" charset=\"utf-8\"></script>\n",
    "\n",
    "Note these are the expected numbers for consecutive rings.\n",
    "\n",
    "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](http://www.4dsolutions.net/ocn/numeracy0.html) and then render in [POV-Ray](http://www.povray.org).  We just want to look at the edges and carbon atom vertexes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sph = \"\"\"sphere { %s 0.1 texture { pigment { color rgb <1,0,0> } } }\"\"\"\n",
    "cyl = \"\"\"cylinder { %s %s 0.05 texture { pigment { color rgb <1.0, 0.65, 0.0> } } }\"\"\"\n",
    "\n",
    "def make_graphene(fname=\"../c6xty/graphene.pov\", append=True):\n",
    "    \"\"\"\n",
    "    Scan through carbons, edges, converting to XYZ and embedding\n",
    "    in POV-Ray Scene Description Language\n",
    "    \"\"\"\n",
    "    if append:\n",
    "        pov = open(fname, \"a\")\n",
    "    else:\n",
    "        pov = open(fname, \"w\")\n",
    "\n",
    "    # graphene will be included as a single object in the \n",
    "    # parent povray script, where lighting, camera position,\n",
    "    # and background are defined\n",
    "    \n",
    "    print(\"#declare graphene = union{\", file=pov)\n",
    "    for atom in carbons:\n",
    "        v = Qvector(atom).xyz()\n",
    "        s = sph % \"<{0.x}, {0.y}, {0.z}>\".format(v)\n",
    "        print(s, file=pov)\n",
    "    for bond in edges:\n",
    "        v0, v1 = bond\n",
    "        v0 = Qvector(v0).xyz()\n",
    "        v1 = Qvector(v1).xyz()\n",
    "        c = cyl % (\"<{0.x}, {0.y}, {0.z}>\".format(v0), \"<{0.x}, {0.y}, {0.z}>\".format(v1))\n",
    "        print(c, file=pov)\n",
    "    print(\"}\\n\", file=pov)\n",
    "        \n",
    "make_graphene(append=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Success!\n",
    "\n",
    "<a data-flickr-embed=\"true\"  href=\"https://www.flickr.com/photos/kirbyurner/44743063302/in/dateposted-public/\" title=\"Another Test\"><img src=\"https://farm2.staticflickr.com/1877/44743063302_e7db33cdea_b.jpg\" width=\"1024\" height=\"768\" alt=\"Another Test\"></a><script async src=\"//embedr.flickr.com/assets/client-code.js\" charset=\"utf-8\"></script>\n",
    "\n",
    "## ADDENDUM\n",
    "\n",
    "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.  \n",
    "\n",
    "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).\n",
    "\n",
    "<a data-flickr-embed=\"true\"  href=\"https://www.flickr.com/photos/kirbyurner/44174173534/in/dateposted-public/\" title=\"Better Sunlight\"><img src=\"https://farm2.staticflickr.com/1958/44174173534_2773c5f9cf_z.jpg\" width=\"640\" height=\"480\" alt=\"Better Sunlight\"></a><script async src=\"//embedr.flickr.com/assets/client-code.js\" charset=\"utf-8\"></script>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}