3DNPs.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b617cc8b-67f4-448d-8416-3ff073e18fba",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/home/romuald/ENSEIGNEMENT/0-JupyterNotebooks/pyPhysChem\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<style>\n",
       "h1, h2, h3, h4, h5, h6 {\n",
       "  font-family: Verdana, \"DejaVu Sans\", \"Bitstream Vera Sans\", Geneva, sans-serif;\n",
       "  font-weight: bold;\n",
       "}\n",
       "body, intro, introE, rq, rqE, ex, exE, app, appE, sol  {\n",
       "  font-family: Verdana, \"DejaVu Sans\", \"Bitstream Vera Sans\", Geneva, sans-serif;\n",
       "  font-weight: 400;\n",
       "  font-size: 12px;\n",
       "}\n",
       "h1 {\n",
       "  border: 0 solid #333;\n",
       "  padding: 30px ;\n",
       "  color: white;\n",
       "  background: #b11d01;\n",
       "  text-align: center;\n",
       "}\n",
       "h2 {\n",
       "  border: 3px solid #333;\n",
       "  padding: 18px ;\n",
       "  color: #b11d01;\n",
       "  background: #ffffff;\n",
       "  text-align: center;\n",
       "}\n",
       "h3 {\n",
       "  border: 0 solid #333;\n",
       "  padding: 12px ;\n",
       "  color: #000000;\n",
       "  background: #c1c1c1;\n",
       "  text-align: left;\n",
       "}\n",
       "h4 {\n",
       "  border: 0 solid #333;\n",
       "  padding: 2px ;\n",
       "  color: #000000;\n",
       "  background: #d9fffc;\n",
       "  text-align: left;\n",
       "}\n",
       "h5 {\n",
       "  border: 1px solid #333;\n",
       "  padding: 2px ;\n",
       "  color: #000000;\n",
       "  background: #ffffff;\n",
       "  text-align: left;\n",
       "}\n",
       ".introT::before {    \n",
       "    content: attr(title);\n",
       "    background-color: #cecece;\n",
       "    margin-left:-0.5em;\n",
       "    margin-right:-0.5em;\n",
       "    padding-left:0.5em;\n",
       "    padding-right:0.5em;\n",
       "    font-weight: 600;\n",
       "    display: block;\n",
       "}\n",
       ".introT {    \n",
       "    background-color: #cecece80;\n",
       "    border-color: #969696;\n",
       "    border-left: 5px solid #969696;\n",
       "    padding: 0.5em;\n",
       "}\n",
       ".intro {    \n",
       "    background-color: #cecece80;\n",
       "    border-color: #969696;\n",
       "    border-left: 5px solid #969696;\n",
       "    padding: 0.5em;\n",
       "}\n",
       ".introE {    \n",
       "    background-color: #cecece80;\n",
       "    border-color: #969696;\n",
       "    border-left: 5px solid #969696;\n",
       "    padding: 0.5em;\n",
       "    color : #117996;\n",
       "}\n",
       ".rq {    \n",
       "    background-color: #fcf2f2;\n",
       "    border-color: #dFb5b4;\n",
       "    border-left: 5px solid #dfb5b4;\n",
       "    padding: 0.5em;\n",
       "}\n",
       ".rqE::before {    \n",
       "    background-color: #fcd3d3;\n",
       "    color: #ff0000;\n",
       "    content:\"Warning\";\n",
       "    margin-left:-0.5em;\n",
       "    margin-right:-0.5em;\n",
       "    padding-left:0.5em;\n",
       "    padding-right:0.5em;\n",
       "    font-weight: 600;\n",
       "    display: block;\n",
       "}\n",
       ".rqE {    \n",
       "    background-color: #fcd3d380;\n",
       "    border-color: #dFb5b4;\n",
       "    border-left: 5px solid #dfb5b4;\n",
       "    padding: 0.5em;\n",
       "}\n",
       ".rqT::before {    \n",
       "    background-color: #fcd3d3;\n",
       "    color: #ff0000;\n",
       "    content: attr(title);\n",
       "    margin-left:-0.5em;\n",
       "    margin-right:-0.5em;\n",
       "    padding-left:0.5em;\n",
       "    padding-right:0.5em;\n",
       "    font-weight: 600;\n",
       "    display: block;\n",
       "}\n",
       ".rqT {    \n",
       "    background-color: #fcd3d380;\n",
       "    border-color: #dFb5b4;\n",
       "    border-left: 5px solid #dfb5b4;\n",
       "    padding: 0.5em;\n",
       "}\n",
       ".exold {    \n",
       "    background-color: #b2dbea80;\n",
       "    border-color: #0055ff;\n",
       "    border-left: 10px solid #0055ff;\n",
       "    padding: 0.5em;\n",
       "}\n",
       ".ex {    \n",
       "    background-color: #b2dbea80;\n",
       "    padding: 0.5em;\n",
       "    padding-top: 0em;\n",
       "    position:relative;\n",
       "}\n",
       ".ex::before {\n",
       "    background-color: #b2dbea;\n",
       "    content:\"Exercice\";\n",
       "    margin-left:-0.5em;\n",
       "    margin-right:-0.5em;\n",
       "    padding-left:0.5em;\n",
       "    padding-right:0.5em;\n",
       "    font-weight: 600;\n",
       "    display: block;\n",
       "}\n",
       ".exE {    \n",
       "    background-color: #b2dbea80;\n",
       "    padding: 0.5em;\n",
       "    padding-top: 0em;\n",
       "    position:relative;\n",
       "    color : #117996;\n",
       "}\n",
       ".exE::before {\n",
       "    background-color: #b2dbea;\n",
       "    content:\"Exercise\";\n",
       "    margin-left:-0.5em;\n",
       "    margin-right:-0.5em;\n",
       "    padding-left:0.5em;\n",
       "    padding-right:0.5em;\n",
       "    font-weight: 600;\n",
       "    display: block;\n",
       "    color : #117996;\n",
       "}\n",
       ".app {    \n",
       "    background-color: #b2dbea80;\n",
       "    padding: 0.5em;\n",
       "    padding-top: 0em;\n",
       "    position:relative;\n",
       "}\n",
       ".app::before {\n",
       "    background-color: #b2dbea;\n",
       "    content:\"Application\";\n",
       "    margin-left:-0.5em;\n",
       "    margin-right:-0.5em;\n",
       "    padding-left:0.5em;\n",
       "    padding-right:0.5em;\n",
       "    font-weight: 600;\n",
       "    display: block;\n",
       "}\n",
       ".appE {    \n",
       "    background-color: #b2dbea80;\n",
       "    padding: 0.5em;\n",
       "    padding-top: 0em;\n",
       "    color : #117996;\n",
       "    position:relative;\n",
       "}\n",
       ".appE::before {\n",
       "    background-color: #b2dbea;\n",
       "    content:\"Application\";\n",
       "    margin-left:-0.5em;\n",
       "    margin-right:-0.5em;\n",
       "    padding-left:0.5em;\n",
       "    padding-right:0.5em;\n",
       "    font-weight: 600;\n",
       "    color : #117996;\n",
       "    display: block;\n",
       "}\n",
       ".sol {    \n",
       "    background-color: #bbeab880;\n",
       "    padding: 0.5em;\n",
       "    padding-top: 0em;\n",
       "    position:relative;\n",
       "}\n",
       ".sol::before {\n",
       "    background-color: #bbeab8;\n",
       "    content:\"Answer\";\n",
       "    margin-left:-0.5em;\n",
       "    margin-right:-0.5em;\n",
       "    padding-left:0.5em;\n",
       "    padding-right:0.5em;\n",
       "    font-weight: 600;\n",
       "    display: block;\n",
       "}\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "**Start at:** Friday 06 October 2023, 15:58:56  \n",
       "**Hostname:** insa-12842 (Linux)"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "<p style=\"text-align: center\"><img width=\"800px\" src=\"./config/svg/pyPhysChemBanner.svg\" style=\"margin-left:auto; margin-right:auto\"/></p>"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import os\n",
    "import sys\n",
    "\n",
    "print(os.getcwd())\n",
    "cwd0 = './config/'\n",
    "sys.path.append(cwd0)\n",
    "\n",
    "import visualID_Eng as vID\n",
    "from visualID_Eng import  fg, hl, bg\n",
    "vID.init(cwd0)\n",
    "import tools4pyPhysChem as t4pPC"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2731bde-6a3e-459b-9e1b-9851717b0888",
   "metadata": {},
   "source": [
    "# 3D NPs. Generation of coordinates"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "127eda72-cbec-44fa-baf4-9b5fbea52900",
   "metadata": {},
   "source": [
    "## Utils"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8895ee46-f0ae-4dec-a6b0-8dfced0d0d55",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from ase.io import write\n",
    "from ase import Atoms\n",
    "from ase.visualize import view"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3db5d65d-6c32-4f86-9b9b-6ff8f46c79e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def RAB(coord,a,b):\n",
    "    import numpy as np\n",
    "    \"\"\"calculate the interatomic distance between two atoms a and b\"\"\"\n",
    "    r = np.sqrt((coord[a][0]-coord[b][0])**2 + (coord[a][1]-coord[b][1])**2 + (coord[a][2]-coord[b][2])**2)\n",
    "    return r\n",
    "\n",
    "def Rbetween2Points(p1,p2):\n",
    "    import numpy as np\n",
    "    r = np.sqrt((p1[0]-p2[0])**2 + (p1[1]-p2[1])**2 + (p1[2]-p2[2])**2)\n",
    "    return r\n",
    "\n",
    "def vector(coord,a,b):\n",
    "    import numpy as np\n",
    "    v = [coord[b][0]-coord[a][0],coord[b][1]-coord[a][1],coord[b][2]-coord[a][2]]\n",
    "    v = np.array(v)\n",
    "    return v\n",
    "\n",
    "def vectorBetween2Points(p1,p2):\n",
    "    import numpy as np\n",
    "    v = [p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]]\n",
    "    v = np.array(v)\n",
    "    return v\n",
    "\n",
    "def coord2xyz(coord):\n",
    "    x = np.array(coord)[:,0]\n",
    "    y = np.array(coord)[:,1]\n",
    "    z = np.array(coord)[:,2]\n",
    "    return x,y,z\n",
    "\n",
    "def vertex(x, y, z, scale):\n",
    "    \"\"\" Return vertex coordinates fixed to the unit sphere \"\"\"\n",
    "    length = np.sqrt(x**2 + y**2 + z**2)\n",
    "    return [(i * scale) / length for i in (x,y,z)]\n",
    "\n",
    "    \n",
    "def RadiusSphereAfterV(V):\n",
    "    return (3*V/(4*np.pi))**(1/3)\n",
    "\n",
    "def MakeFaceCoord(Rnn,f,coord,nAtomsOnFaces,coordFaceAt):\n",
    "    # adapt these two lines for square faces or other\n",
    "    # the idea here is to interpolate between edge atoms of two relevant edges\n",
    "    # (for example two opposite edges of a squared face)\n",
    "    # be careful of the vectors orientation of the edges!\n",
    "    if (len(f) == 3):  #triangular facet\n",
    "        edge1 = [f[1],f[0]]\n",
    "        edge2 = [f[1],f[2]]\n",
    "    if (len(f) == 5):  #pentagonal facet\n",
    "        edge1 = [f[1],f[0]]\n",
    "        edge2 = [f[1],f[2]]\n",
    "    nAtomsOnEdges = int((RAB(coord,f[1],f[0])+1e-6)/Rnn) - 1\n",
    "    nIntervalsE = nAtomsOnEdges + 1\n",
    "    #print(\"shell = \",i,\"nAtomsOnEdges = \",nAtomsOnEdges)\n",
    "    for n in range(nAtomsOnEdges):\n",
    "        CoordAtomOnEdge1 = coord[edge1[0]]+vector(coord,edge1[0],edge1[1])*(n+1) / nIntervalsE\n",
    "        CoordAtomOnEdge2 = coord[edge2[0]]+vector(coord,edge2[0],edge2[1])*(n+1) / nIntervalsE\n",
    "        #print(CoordAtomOnEdge1)\n",
    "        #print(CoordAtomOnEdge2)\n",
    "        distBetween2EdgeAtoms = Rbetween2Points(CoordAtomOnEdge1,CoordAtomOnEdge2)\n",
    "        #print(\"R = \",distBetween2EdgeAtoms)\n",
    "        nAtomsBetweenEdges = int((distBetween2EdgeAtoms+1e-6)/Rnn) - 1\n",
    "        nIntervalsF = nAtomsBetweenEdges + 1\n",
    "        for m in range(nAtomsBetweenEdges):\n",
    "            coordFaceAt.append(CoordAtomOnEdge1 + vectorBetween2Points(CoordAtomOnEdge1,CoordAtomOnEdge2)*(m+1) / nIntervalsF)\n",
    "            nAtomsOnFaces += 1\n",
    "    return nAtomsOnFaces,coordFaceAt\n",
    "\n",
    "def centerOfGravity(c,select):\n",
    "    nselect = len(select)\n",
    "    xg = 0\n",
    "    yg = 0\n",
    "    zg = 0\n",
    "    for at in select:\n",
    "        xg += c[at][0]\n",
    "        yg += c[at][1]\n",
    "        zg += c[at][2]\n",
    "    cog = [xg/nselect, yg/nselect, zg/nselect]\n",
    "    return cog\n",
    "\n",
    "# def principal_moment_of_inertia(mol):\n",
    "#     from ase import Atoms\n",
    "#     com = mol.get_center_of_mass()\n",
    "#     m = mol.get_masses()\n",
    "#     nAtoms = len(mol)\n",
    "#     x,y,z = coord2xyz(mol.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdf9253b-5a5f-4c96-9b96-255414d14dbb",
   "metadata": {},
   "source": [
    "## Icosahedra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc3b7ad4-246e-4404-832d-86ce0cc8ef94",
   "metadata": {},
   "source": [
    "### Regular icosahedra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "93c67994-9424-44f3-96fe-6dcc8fac8f4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.cluster import Icosahedron\n",
    "\n",
    "atoms = Icosahedron(\"Au\", 8, 4.1536)\n",
    "atoms.center(vacuum=10)\n",
    "\n",
    "write(\"3DNP-coords/ico.xyz\", atoms)\n",
    "write(\"3DNP-coords/ico.cif\", atoms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "40c85e33-1275-473c-9f8d-f5f475b0380f",
   "metadata": {},
   "outputs": [],
   "source": [
    "class regIC:\n",
    "    import numpy as np\n",
    "    nFaces = 20\n",
    "    nEdges = 30\n",
    "    nVertices = 12\n",
    "    phi = (1 + np.sqrt(5))/2 # golden ratio\n",
    "    edgeLengthF = 1\n",
    "    radiusCSF = np.sqrt(10 + 2*np.sqrt(5))/4\n",
    "    interShellF = 1/radiusCSF\n",
    "#    interShellF = np.sqrt(2*(1-1/np.sqrt(5)))\n",
    "    radiusISF = np.sqrt(3) * (3 + np.sqrt(5))/12\n",
    "  \n",
    "    def __init__(self, Rnn, nShell):\n",
    "        self.Rnn = Rnn\n",
    "        self.nShell = nShell\n",
    "        self.nAtoms = 0\n",
    "        self.nAtomsPerShell = [0]\n",
    "        self.interShellDistance = self.Rnn / self.interShellF\n",
    "          \n",
    "    def __str__(self):\n",
    "        return(f\"Regular icosahedron with {self.nShell} shell(s) and Rnn = {self.Rnn}\")\n",
    "    \n",
    "    def nAtomsF(self,i):\n",
    "        return (10*i**3 + 11*i)//3 + 5*i**2 + 1\n",
    "    \n",
    "    def nAtomsPerShellAnalytic(self):\n",
    "        n = []\n",
    "        Sum = 0\n",
    "        for i in range(self.nShell+1):\n",
    "            Sum = sum(n)\n",
    "            ni = self.nAtomsF(i)\n",
    "            n.append(ni-Sum)\n",
    "        return n\n",
    "    \n",
    "    def nAtomsAnalytic(self):\n",
    "        n = self.nAtomsF(self.nShell)\n",
    "        return n\n",
    "    \n",
    "    def edgeLength(self):\n",
    "        return self.Rnn*self.nShell\n",
    "\n",
    "    def radiusCircumscribedSphere(self):\n",
    "        return self.radiusCSF*self.edgeLength()\n",
    "\n",
    "    def radiusInscribedSphere(self):\n",
    "        return self.radiusISF*self.edgeLength()\n",
    "\n",
    "    def area(self):\n",
    "        el = self.edgeLength()\n",
    "        return 5 * el**2 * np.sqrt(3)\n",
    "    \n",
    "    def volume(self):\n",
    "        el = self.edgeLength()\n",
    "        return 5 * el**3 * (3 + np.sqrt(5))/12\n",
    "    \n",
    "    def MakeVertices(self,i):\n",
    "        import sys\n",
    "        \"\"\"\n",
    "        input:\n",
    "            - i = index of the shell\n",
    "        returns:\n",
    "            - CoordVertices = the 12 vertex coordinates of the ith shell of an icosahedron\n",
    "            - edges = indexes of the 30 edges\n",
    "            - faces = indexes of the 20 faces \n",
    "        \"\"\"\n",
    "        if (i == 0):\n",
    "            CoordVertices = [0., 0., 0.]\n",
    "            edges = []\n",
    "            faces = []\n",
    "        elif (i > self.nShell):\n",
    "            sys.exit(f\"icoreg.MakeVertices(i) is called with i = {i} > nShell = {self.nShell}\")\n",
    "        else:\n",
    "            phi = self.phi\n",
    "            scale = self.interShellDistance * i\n",
    "            CoordVertices = [ vertex(-1, phi, 0, scale),\\\n",
    "                       vertex( 1, phi, 0, scale),\\\n",
    "                       vertex(-1, -phi, 0, scale),\\\n",
    "                       vertex( 1, -phi, 0, scale),\\\n",
    "                       vertex(0, -1, phi, scale),\\\n",
    "                       vertex(0, 1, phi, scale),\\\n",
    "                       vertex(0, -1, -phi, scale),\\\n",
    "                       vertex(0, 1, -phi, scale),\\\n",
    "                       vertex( phi, 0, -1, scale),\\\n",
    "                       vertex( phi, 0, 1, scale),\\\n",
    "                       vertex(-phi, 0, -1, scale),\\\n",
    "                       vertex(-phi, 0, 1, scale) ]\n",
    "            edges = [( 1, 0), ( 3, 2), ( 4, 2), ( 4, 3), ( 5, 0), ( 5, 1), ( 5, 4), ( 6, 2), ( 6, 3), ( 7, 0),\\\n",
    "                  ( 7, 1), ( 7, 6), ( 8, 1), ( 8, 3), ( 8, 6), ( 8, 7), ( 9, 1), ( 9, 3), ( 9, 4), ( 9, 5),\\\n",
    "                  ( 9, 8), (10, 0), (10, 2), (10, 6), (10, 7), (11, 0), (11, 2), (11, 4), (11, 5), (11,10),]\n",
    "            faces = [(7,0,1),(7,1,8),(7,8,6),(7,6,10),(7,10,0),\\\n",
    "                     (0,11,5),(0,5,1),(1,5,9),(1,8,9),(8,9,3),(8,3,6),(6,3,2),(6,10,2),(10,2,11),(10,0,11),\\\n",
    "                     (4,2,3),(4,3,9),(4,9,5),(4,5,11),(4,11,2)]\n",
    "            edges = np.array(edges)\n",
    "            CoordVertices = np.array(CoordVertices)\n",
    "            faces = np.array(faces)\n",
    "        return CoordVertices, edges, faces\n",
    "\n",
    "    def coords(self):\n",
    "        # central atom = \"1st shell\"\n",
    "        c = [[0., 0., 0.]]\n",
    "        self.nAtoms = 1\n",
    "        self.nAtomsPerShell = [0]\n",
    "        self.nAtomsPerShell[0] = 1\n",
    "        indexVertexAtoms = []\n",
    "        indexEdgeAtoms = []\n",
    "        indexFaceAtoms = []\n",
    "        for i in range(1,self.nShell+1):\n",
    "            # vertices\n",
    "            nAtoms0 = self.nAtoms\n",
    "            cshell, E, F = self.MakeVertices(i)\n",
    "            self.nAtoms += self.nVertices\n",
    "            self.nAtomsPerShell.append(self.nVertices)\n",
    "            c.extend(cshell.tolist())\n",
    "            indexVertexAtoms.extend(range(nAtoms0,self.nAtoms))\n",
    "\n",
    "            # intermediate atoms on edges e\n",
    "            nAtoms0 = self.nAtoms\n",
    "            Rvv = RAB(cshell,E[0,0],E[0,1]) #distance between two vertex atoms\n",
    "            nAtomsOnEdges = int((Rvv+1e-6) / self.Rnn)-1\n",
    "            nIntervals = nAtomsOnEdges + 1\n",
    "            #print(\"nAtomsOnEdges = \",nAtomsOnEdges)\n",
    "            coordEdgeAt = []\n",
    "            for n in range(nAtomsOnEdges):\n",
    "                for e in E:\n",
    "                    a = e[0]\n",
    "                    b = e[1]\n",
    "                    coordEdgeAt.append(cshell[a]+vector(cshell,a,b)*(n+1) / nIntervals)\n",
    "            self.nAtomsPerShell[i] += nAtomsOnEdges * len(E) # number of edges x nAtomsOnEdges\n",
    "            self.nAtoms += nAtomsOnEdges * len(E)\n",
    "            c.extend(coordEdgeAt)\n",
    "            indexEdgeAtoms.extend(range(nAtoms0,self.nAtoms))\n",
    "            \n",
    "            # now, facet atoms\n",
    "            coordFaceAt = []\n",
    "            nAtomsOnFaces = 0\n",
    "            nAtoms0 = self.nAtoms\n",
    "            for f in F:\n",
    "                nAtomsOnFaces,coordFaceAt = MakeFaceCoord(self.Rnn,f,cshell,nAtomsOnFaces,coordFaceAt)\n",
    "            self.nAtomsPerShell[i] += nAtomsOnFaces\n",
    "            self.nAtoms += nAtomsOnFaces\n",
    "            c.extend(coordFaceAt)\n",
    "            indexFaceAtoms.extend(range(nAtoms0,self.nAtoms))\n",
    "            \n",
    "        # print(indexVertexAtoms)\n",
    "        # print(indexEdgeAtoms)\n",
    "        # print(indexFaceAtoms)\n",
    "        return c,[indexVertexAtoms,indexEdgeAtoms,indexFaceAtoms]\n",
    "    \n",
    "    def prop(self):\n",
    "        print(self)\n",
    "        print(\"number of vertices = \",self.nVertices)\n",
    "        print(\"number of edges = \",self.nEdges)\n",
    "        print(\"number of faces = \",self.nFaces)\n",
    "        print(\"phi = \",self.phi)\n",
    "        print(f\"intershell factor = {self.interShellF:.2f}\")\n",
    "        print(f\"nearest neighbour distance = {self.Rnn:.2f} Å\")\n",
    "        print(f\"intershell distance = {self.interShellDistance:.2f} Å\")\n",
    "        print(f\"edge length = {self.edgeLength()*0.1:.2f} nm\")\n",
    "        print(f\"radius after volume = {RadiusSphereAfterV(self.volume()*1e-3):.2f} nm\")\n",
    "        print(f\"radius of the circumscribed sphere = {self.radiusCircumscribedSphere()*0.1:.2f} nm\")\n",
    "        print(f\"radius of the inscribed sphere = {self.radiusInscribedSphere()*0.1:.2f} nm\")\n",
    "        print(f\"area = {self.area()*1e-2:.1f} nm2\")\n",
    "        print(f\"volume = {self.volume()*1e-3:.1f} nm3\")\n",
    "        print(\"number of atoms per shell = \",self.nAtomsPerShellAnalytic())\n",
    "        print(\"total number of atoms = \",self.nAtomsAnalytic())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "975eeda4-d8d3-42fb-bed6-6350e05fcb1b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Regular icosahedron with 4 shell(s) and Rnn = 2.7\n",
      "number of vertices =  12\n",
      "number of edges =  30\n",
      "number of faces =  20\n",
      "phi =  1.618033988749895\n",
      "intershell factor = 1.05\n",
      "nearest neighbour distance = 2.70 Å\n",
      "intershell distance = 2.57 Å\n",
      "edge length = 1.08 nm\n",
      "radius after volume = 0.87 nm\n",
      "radius of the circumscribed sphere = 1.03 nm\n",
      "radius of the inscribed sphere = 0.82 nm\n",
      "area = 10.1 nm2\n",
      "volume = 2.7 nm3\n",
      "number of atoms per shell =  [1, 12, 42, 92, 162]\n",
      "total number of atoms =  309\n",
      "309\n",
      "[1, 12, 42, 92, 162]\n",
      "Moments of inertia = 2341672.64 2341672.64 2341672.64 amu.Å2\n",
      "Moments of inertia / M = 38.47 38.47 38.47 amu.Å2\n",
      "Size of the ellipsoid = 2.77 2.77 2.77 nm\n",
      "\n",
      "Regular icosahedron with 10 shell(s) and Rnn = 2.7\n",
      "number of vertices =  12\n",
      "number of edges =  30\n",
      "number of faces =  20\n",
      "phi =  1.618033988749895\n",
      "intershell factor = 1.05\n",
      "nearest neighbour distance = 2.70 Å\n",
      "intershell distance = 2.57 Å\n",
      "edge length = 2.70 nm\n",
      "radius after volume = 2.17 nm\n",
      "radius of the circumscribed sphere = 2.57 nm\n",
      "radius of the inscribed sphere = 2.04 nm\n",
      "area = 63.1 nm2\n",
      "volume = 42.9 nm3\n",
      "number of atoms per shell =  [1, 12, 42, 92, 162, 252, 362, 492, 642, 812, 1002]\n",
      "total number of atoms =  3871\n",
      "3871\n",
      "[1, 12, 42, 92, 162, 252, 362, 492, 642, 812, 1002]\n",
      "Moments of inertia / M = 210.25 210.25 210.25 amu.Å2\n",
      "Size of the ellipsoid = 6.48 6.48 6.48 nm\n",
      "\n",
      "Regular icosahedron with 20 shell(s) and Rnn = 2.7\n",
      "number of vertices =  12\n",
      "number of edges =  30\n",
      "number of faces =  20\n",
      "phi =  1.618033988749895\n",
      "intershell factor = 1.05\n",
      "nearest neighbour distance = 2.70 Å\n",
      "intershell distance = 2.57 Å\n",
      "edge length = 5.40 nm\n",
      "radius after volume = 4.34 nm\n",
      "radius of the circumscribed sphere = 5.14 nm\n",
      "radius of the inscribed sphere = 4.08 nm\n",
      "area = 252.5 nm2\n",
      "volume = 343.5 nm3\n",
      "number of atoms per shell =  [1, 12, 42, 92, 162, 252, 362, 492, 642, 812, 1002, 1212, 1442, 1692, 1962, 2252, 2562, 2892, 3242, 3612, 4002]\n",
      "total number of atoms =  28741\n",
      "28741\n",
      "[1, 12, 42, 92, 162, 252, 362, 492, 642, 812, 1002, 1212, 1442, 1692, 1962, 2252, 2562, 2892, 3242, 3612, 4002]\n"
     ]
    }
   ],
   "source": [
    "IC1 = regIC(2.7,4)\n",
    "IC1.prop()\n",
    "c1,i1 = IC1.coords()\n",
    "print(IC1.nAtoms)\n",
    "print(IC1.nAtomsPerShell)\n",
    "ico1 = Atoms('Au'*IC1.nAtoms, positions=c1)\n",
    "write(\"3DNP-coords/ico1.xyz\", ico1)\n",
    "ico1.moi = ico1.get_moments_of_inertia() # in amu*angstrom**2\n",
    "print(f\"Moments of inertia = {ico1.moi[0]:.2f} {ico1.moi[1]:.2f} {ico1.moi[2]:.2f} amu.Å2\")\n",
    "ico1.masses = ico1.get_masses()\n",
    "ico1.M = ico1.masses.sum()\n",
    "ico1.moiM = ico1.moi/ico1.M\n",
    "print(f\"Moments of inertia / M = {ico1.moiM[0]:.2f} {ico1.moiM[1]:.2f} {ico1.moiM[2]:.2f} amu.Å2\")\n",
    "ico1.dim = 2*np.sqrt(5*ico1.moiM)\n",
    "print(f\"Size of the ellipsoid = {ico1.dim[0]*0.1:.2f} {ico1.dim[1]*0.1:.2f} {ico1.dim[2]*0.1:.2f} nm\")\n",
    "\n",
    "#view(ico1)\n",
    "write(\"3DNP-coords/ico1.xyz\", ico1)\n",
    "print()\n",
    "IC2 = regIC(2.7,10)\n",
    "IC2.prop()\n",
    "c2,i2 = IC2.coords()\n",
    "print(IC2.nAtoms)\n",
    "print(IC2.nAtomsPerShell)\n",
    "ico2 = Atoms('Au'*IC2.nAtoms, positions=c2)\n",
    "#view(ico2)\n",
    "write(\"3DNP-coords/ico2.xyz\", ico2)\n",
    "ico2.moi = ico2.get_moments_of_inertia() # in amu*angstrom**2\n",
    "ico2.masses = ico2.get_masses()\n",
    "ico2.M = ico2.masses.sum()\n",
    "ico2.moiM = ico2.moi/ico2.M\n",
    "print(f\"Moments of inertia / M = {ico2.moiM[0]:.2f} {ico2.moiM[1]:.2f} {ico2.moiM[2]:.2f} amu.Å2\")\n",
    "ico2.dim = 2*np.sqrt(5*ico2.moiM)\n",
    "print(f\"Size of the ellipsoid = {ico2.dim[0]*0.1:.2f} {ico2.dim[1]*0.1:.2f} {ico2.dim[2]*0.1:.2f} nm\")\n",
    "\n",
    "print()\n",
    "IC3 = regIC(2.7,20)\n",
    "IC3.prop()\n",
    "c3,i3 = IC3.coords()\n",
    "print(IC3.nAtoms)\n",
    "print(IC3.nAtomsPerShell)\n",
    "ico3 = Atoms('Au'*IC3.nAtoms, positions=c3)\n",
    "#view(ico2)\n",
    "write(\"3DNP-coords/ico3.xyz\", ico3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f704abb-8484-4363-995c-6b858f856ff4",
   "metadata": {},
   "source": [
    "## Tetrahedra"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8a553c0-47d2-4d5e-ae1f-a9e9a44536e6",
   "metadata": {},
   "source": [
    "### Regular tetrahedron"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d95b31bf-7738-4808-90c3-876880ddfff2",
   "metadata": {},
   "outputs": [],
   "source": [
    "class regTd:\n",
    "    import numpy as np\n",
    "    nFaces = 4\n",
    "    nEdges = 6\n",
    "    nVertices = 4\n",
    "    edgeLengthF = 2\n",
    "    heightOfPyramid = edgeLengthF * np.sqrt(2/3)\n",
    "    radiusCSF = edgeLengthF * np.sqrt(3/8)\n",
    "    radiusISF = edgeLengthF/np.sqrt(24)\n",
    "    radiusMSF = edgeLengthF/np.sqrt(8)\n",
    "    fveAngle = np.rad2deg(np.arccos(1/np.sqrt(3)))\n",
    "    fefAngle = np.rad2deg(np.arccos(1/3))\n",
    "    vcvAngle = np.rad2deg(np.arccos(-1/3))\n",
    "  \n",
    "    def __init__(self, Rnn, nSlice):\n",
    "        self.Rnn = Rnn\n",
    "        self.nSlice = nSlice\n",
    "        self.nAtoms = 0\n",
    "        self.nAtomsPerSlice = []\n",
    "        self.interSliceDistance = self.Rnn / self.heightOfPyramid\n",
    "          \n",
    "    def __str__(self):\n",
    "        return(f\"Regular dodecahedron with {self.nSlice} slice(s) and Rnn = {self.Rnn}\")\n",
    "    \n",
    "    def nAtomsF(self,i):\n",
    "        return i**3/6 + i**2 + 11*i/6 + 1\n",
    "    \n",
    "    def nAtomsPerSliceAnalytic(self):\n",
    "        n = []\n",
    "        Sum = 0\n",
    "        for i in range(self.nSlice+1):\n",
    "            Sum = sum(n)\n",
    "            ni = self.nAtomsF(i)\n",
    "            n.append(ni-Sum)\n",
    "            print(i,ni,Sum,n)\n",
    "        return n\n",
    "    \n",
    "    def nAtomsAnalytic(self):\n",
    "        n = self.nAtomsF(self.nSlice)\n",
    "        return n\n",
    "    \n",
    "    def edgeLength(self):\n",
    "        return self.Rnn*self.nSlice\n",
    "\n",
    "    def radiusCircumscribedSphere(self):\n",
    "        return self.radiusCSF*self.edgeLength()\n",
    "\n",
    "    def radiusInscribedSphere(self):\n",
    "        return self.radiusISF*self.edgeLength()\n",
    "\n",
    "    def radiusMidSphere(self):\n",
    "        return self.radiusMSF*self.edgeLength()\n",
    "\n",
    "    def area(self):\n",
    "        el = self.edgeLength()\n",
    "        return el**2*np.sqrt(3)\n",
    "    \n",
    "    def volume(self):\n",
    "        el = self.edgeLength()\n",
    "        return el**3/(6*np.sqrt(2)) \n",
    "\n",
    "    def MakeVertices(self,i):\n",
    "        import sys\n",
    "        \"\"\"\n",
    "        input:\n",
    "            - i = index of the slice\n",
    "        returns:\n",
    "            - CoordVertices = the 4 vertex coordinates of the ith slice of a tetrahedron\n",
    "            - edges = indexes of the 6 edges\n",
    "            - faces = indexes of the 4 faces \n",
    "        \"\"\"\n",
    "        if (i > self.nSlice):\n",
    "            sys.exit(f\"icoreg.MakeVertices(i) is called with i = {i} > nSlice = {self.nSlice}\")\n",
    "        else:\n",
    "            scale = self.interSliceDistance * i\n",
    "            CoordVertices = [vertex(1, 1, 1, scale),\\\n",
    "                             vertex(1, -1, -1, scale),\\\n",
    "                             vertex(-1, 1, -1, scale),\\\n",
    "                             vertex(-1, -1, 1, scale)]\n",
    "            edges = [(0,1), (0,2), (0,3), (1,2), (1,3), (2,3)]\n",
    "            faces = [(0,2,1),(0,1,3),(0,3,2),(1,2,3)]\n",
    "            edges = np.array(edges)\n",
    "            CoordVertices = np.array(CoordVertices)\n",
    "            faces = np.array(faces)\n",
    "        return CoordVertices, edges, faces\n",
    "\n",
    "    def coords(self):\n",
    "        c = []\n",
    "        print(self.nAtomsPerSlice)\n",
    "        indexVertexAtoms = []\n",
    "        indexEdgeAtoms = []\n",
    "        indexFaceAtoms = []\n",
    "\n",
    "        # vertices\n",
    "        nAtoms0 = 0\n",
    "        self.nAtoms += self.nVertices\n",
    "        cSlice, E, F = self.MakeVertices(self.nSlice)\n",
    "        c.extend(cSlice.tolist())\n",
    "        indexVertexAtoms.extend(range(nAtoms0,self.nAtoms))\n",
    "\n",
    "        # intermediate atoms on edges e\n",
    "        nAtoms0 = self.nAtoms\n",
    "        Rvv = RAB(cSlice,E[0,0],E[0,1]) #distance between two vertex atoms\n",
    "        nAtomsOnEdges = int((Rvv+1e-6) / self.Rnn)-1\n",
    "        nIntervals = nAtomsOnEdges + 1\n",
    "        #print(\"nAtomsOnEdges = \",nAtomsOnEdges)\n",
    "        coordEdgeAt = []\n",
    "        for n in range(nAtomsOnEdges):\n",
    "            for e in E:\n",
    "                a = e[0]\n",
    "                b = e[1]\n",
    "                coordEdgeAt.append(cSlice[a]+vector(cSlice,a,b)*(n+1) / nIntervals)\n",
    "        self.nAtoms += nAtomsOnEdges * len(E)\n",
    "        c.extend(coordEdgeAt)\n",
    "        indexEdgeAtoms.extend(range(nAtoms0,self.nAtoms))\n",
    "        \n",
    "        # now, facet atoms\n",
    "        coordFaceAt = []\n",
    "        nAtomsOnFaces = 0\n",
    "        nAtoms0 = self.nAtoms\n",
    "        for f in F:\n",
    "            nAtomsOnFaces,coordFaceAt = MakeFaceCoord(self.Rnn,f,c,nAtomsOnFaces,coordFaceAt)\n",
    "        self.nAtoms += nAtomsOnFaces\n",
    "        c.extend(coordFaceAt)\n",
    "        indexFaceAtoms.extend(range(nAtoms0,self.nAtoms))\n",
    "\n",
    "#         for i in range(0,):\n",
    "#             print(i,\"Slice\")\n",
    "#             # vertices\n",
    "#             self.nAtoms += self.nVertices\n",
    "#             self.nAtomsPerSlice.append(self.nVertices)\n",
    "#             c.extend(cSlice.tolist())\n",
    "#             indexVertexAtoms.extend(range(nAtoms0,self.nAtoms))\n",
    "\n",
    "#             # intermediate atoms on edges e\n",
    "#             nAtoms0 = self.nAtoms\n",
    "#             Rvv = RAB(cSlice,E[0,0],E[0,1]) #distance between two vertex atoms\n",
    "#             nAtomsOnEdges = int((Rvv+1e-6) / self.Rnn)-1\n",
    "#             nIntervals = nAtomsOnEdges + 1\n",
    "#             #print(\"nAtomsOnEdges = \",nAtomsOnEdges)\n",
    "#             coordEdgeAt = []\n",
    "#             for n in range(nAtomsOnEdges):\n",
    "#                 for e in E:\n",
    "#                     a = e[0]\n",
    "#                     b = e[1]\n",
    "#                     coordEdgeAt.append(cSlice[a]+vector(cSlice,a,b)*(n+1) / nIntervals)\n",
    "#             self.nAtomsPerSlice[i] += nAtomsOnEdges * len(E) # number of edges x nAtomsOnEdges\n",
    "#             self.nAtoms += nAtomsOnEdges * len(E)\n",
    "#             c.extend(coordEdgeAt)\n",
    "#             indexEdgeAtoms.extend(range(nAtoms0,self.nAtoms))\n",
    "            \n",
    "        return c,[indexVertexAtoms,indexEdgeAtoms,indexFaceAtoms]\n",
    "    \n",
    "    def prop(self):\n",
    "        print(self)\n",
    "        print(\"number of vertices = \",self.nVertices)\n",
    "        print(\"number of edges = \",self.nEdges)\n",
    "        print(\"number of faces = \",self.nFaces)\n",
    "        print(f\"nearest neighbour distance = {self.Rnn:.2f} Å\")\n",
    "        print(f\"edge length = {self.edgeLength()*0.1:.2f} nm\")\n",
    "        print(f\"height of pyramid = {self.heightOfPyramid*0.1:.2f} nm\")\n",
    "        print(f\"radius after volume = {RadiusSphereAfterV(self.volume()*1e-3):.2f} nm\")\n",
    "        print(f\"radius of the circumscribed sphere = {self.radiusCircumscribedSphere()*0.1:.2f} nm\")\n",
    "        print(f\"radius of the inscribed sphere = {self.radiusInscribedSphere()*0.1:.2f} nm\")\n",
    "        print(f\"radius of the midsphere that is tangent to edges = {self.radiusMidSphere()*0.1:.2f} nm\")\n",
    "        print(f\"area = {self.area()*1e-2:.1f} nm2\")\n",
    "        print(f\"volume = {self.volume()*1e-3:.1f} nm3\")\n",
    "        print(f\"face-vertex-edge angle = {self.fveAngle:.1f}°\")\n",
    "        print(f\"face-edge-face (dihedral) angle = {self.fefAngle:.1f}°\")\n",
    "        print(f\"vertex-center-vertex (tetrahedral bond) angle = {self.vcvAngle:.1f}°\")\n",
    "        print(\"number of atoms per slice = \",self.nAtomsPerSliceAnalytic())\n",
    "        print(\"total number of atoms = \",self.nAtomsAnalytic())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "267e77d7-ca75-4838-8244-fb3c543a21c5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Regular dodecahedron with 6 slice(s) and Rnn = 2.7\n",
      "number of vertices =  4\n",
      "number of edges =  6\n",
      "number of faces =  4\n",
      "nearest neighbour distance = 2.70 Å\n",
      "edge length = 1.62 nm\n",
      "height of pyramid = 0.16 nm\n",
      "radius after volume = 0.49 nm\n",
      "radius of the circumscribed sphere = 1.98 nm\n",
      "radius of the inscribed sphere = 0.66 nm\n",
      "radius of the midsphere that is tangent to edges = 1.15 nm\n",
      "area = 4.5 nm2\n",
      "volume = 0.5 nm3\n",
      "face-vertex-edge angle = 54.7°\n",
      "face-edge-face (dihedral) angle = 70.5°\n",
      "vertex-center-vertex (tetrahedral bond) angle = 109.5°\n",
      "0 1.0 0 [1.0]\n",
      "1 4.0 1.0 [1.0, 3.0]\n",
      "2 10.0 4.0 [1.0, 3.0, 6.0]\n",
      "3 20.0 10.0 [1.0, 3.0, 6.0, 10.0]\n",
      "4 35.0 20.0 [1.0, 3.0, 6.0, 10.0, 15.0]\n",
      "5 55.99999999999999 35.0 [1.0, 3.0, 6.0, 10.0, 15.0, 20.999999999999993]\n",
      "6 84.0 55.99999999999999 [1.0, 3.0, 6.0, 10.0, 15.0, 20.999999999999993, 28.000000000000007]\n",
      "number of atoms per slice =  [1.0, 3.0, 6.0, 10.0, 15.0, 20.999999999999993, 28.000000000000007]\n",
      "total number of atoms =  84.0\n",
      "[]\n",
      "Total number of atoms = 74\n",
      "Number of atoms per slice = []\n",
      "Moments of inertia = 353228.03 353228.03 353228.03 amu.Å2\n",
      "Moments of inertia / M = 24.23 24.23 24.23 amu.Å2\n",
      "Size of the ellipsoid = 2.20 2.20 2.20 nm\n"
     ]
    }
   ],
   "source": [
    "Td1 = regTd(2.7,6)\n",
    "Td1.prop()\n",
    "c1,i1 = Td1.coords()\n",
    "print(f\"Total number of atoms = {Td1.nAtoms:.0f}\")\n",
    "print(f\"Number of atoms per slice = {Td1.nAtomsPerSlice}\")\n",
    "Td1 = Atoms('Au'*Td1.nAtoms, positions=c1)\n",
    "write(\"3DNP-coords/Td1.xyz\", Td1)\n",
    "Td1.moi = Td1.get_moments_of_inertia() # in amu*angstrom**2\n",
    "print(f\"Moments of inertia = {Td1.moi[0]:.2f} {Td1.moi[1]:.2f} {Td1.moi[2]:.2f} amu.Å2\")\n",
    "Td1.masses = Td1.get_masses()\n",
    "Td1.M = Td1.masses.sum()\n",
    "Td1.moiM = Td1.moi/Td1.M\n",
    "print(f\"Moments of inertia / M = {Td1.moiM[0]:.2f} {Td1.moiM[1]:.2f} {Td1.moiM[2]:.2f} amu.Å2\")\n",
    "Td1.dim = 2*np.sqrt(5*Td1.moiM)\n",
    "print(f\"Size of the ellipsoid = {Td1.dim[0]*0.1:.2f} {Td1.dim[1]*0.1:.2f} {Td1.dim[2]*0.1:.2f} nm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96801911-7c09-46a8-add7-a8a3ce2d77c9",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Dodecahedra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "3952025c-4561-48dd-b450-3b716cefd0de",
   "metadata": {},
   "outputs": [],
   "source": [
    "class regDD:\n",
    "    import numpy as np\n",
    "    nFaces = 12\n",
    "    nEdges = 30\n",
    "    nVertices = 20\n",
    "    phi = (1 + np.sqrt(5))/2 # golden ratio\n",
    "    edgeLengthF = 1\n",
    "    radiusCSF = np.sqrt(3) * (1 + np.sqrt(5))/4\n",
    "    interShellF = 1/radiusCSF\n",
    "    radiusISF = np.sqrt((5/2) + (11/10)*np.sqrt(5))/2\n",
    "  \n",
    "    def __init__(self, Rnn, nShell):\n",
    "        self.Rnn = Rnn\n",
    "        self.nShell = nShell\n",
    "        self.nAtoms = 0\n",
    "        self.nAtomsPerShell = [0]\n",
    "        self.interShellDistance = self.Rnn / self.interShellF\n",
    "          \n",
    "    def __str__(self):\n",
    "        return(f\"Regular dodecahedron with {self.nShell} shell(s) and Rnn = {self.Rnn}\")\n",
    "    \n",
    "    def nAtomsF(self,i):\n",
    "        return 10*i**3 + 15*i**2 + 7*i + 1\n",
    "    \n",
    "    def nAtomsPerShellAnalytic(self):\n",
    "        n = []\n",
    "        Sum = 0\n",
    "        for i in range(self.nShell+1):\n",
    "            Sum = sum(n)\n",
    "            ni = self.nAtomsF(i)\n",
    "            n.append(ni-Sum)\n",
    "        return n\n",
    "    \n",
    "    def nAtomsAnalytic(self):\n",
    "        n = self.nAtomsF(self.nShell)\n",
    "        return n\n",
    "    \n",
    "    def edgeLength(self):\n",
    "        return self.Rnn*self.nShell\n",
    "\n",
    "    def radiusCircumscribedSphere(self):\n",
    "        return self.radiusCSF*self.edgeLength()\n",
    "\n",
    "    def radiusInscribedSphere(self):\n",
    "        return self.radiusISF*self.edgeLength()\n",
    "\n",
    "    def area(self):\n",
    "        el = self.edgeLength()\n",
    "        return 3 * el**2 * np.sqrt(25 + 10*np.sqrt(5))\n",
    "    \n",
    "    def volume(self):\n",
    "        el = self.edgeLength()\n",
    "        return (15 + 7*np.sqrt(5)) * el**2/4 \n",
    "\n",
    "    def MakeVertices(self,i):\n",
    "        import sys\n",
    "        \"\"\"\n",
    "        input:\n",
    "            - i = index of the shell\n",
    "        returns:\n",
    "            - CoordVertices = the 20 vertex coordinates of the ith shell of a dodecahedron\n",
    "            - edges = indexes of the 30 edges\n",
    "            - faces = indexes of the 12 faces \n",
    "        \"\"\"\n",
    "        if (i == 0):\n",
    "            CoordVertices = [0., 0., 0.]\n",
    "            edges = []\n",
    "            faces = []\n",
    "        elif (i > self.nShell):\n",
    "            sys.exit(f\"icoreg.MakeVertices(i) is called with i = {i} > nShell = {self.nShell}\")\n",
    "        else:\n",
    "            phi = self.phi\n",
    "            scale = self.interShellDistance * i\n",
    "            CoordVertices = [vertex(1, 1, 1, scale),\\\n",
    "                             vertex(-1, 1, 1, scale),\\\n",
    "                             vertex(1, -1, 1, scale),\\\n",
    "                             vertex(1, 1, -1, scale),\\\n",
    "                             vertex(-1, -1, 1, scale),\\\n",
    "                             vertex(-1, 1, -1, scale),\\\n",
    "                             vertex(1, -1, -1, scale),\\\n",
    "                             vertex(-1, -1, -1, scale),\\\n",
    "                             vertex(0, phi, 1/phi, scale),\\\n",
    "                             vertex(0, -phi, 1/phi, scale),\\\n",
    "                             vertex(0, phi, -1/phi, scale),\\\n",
    "                             vertex(0, -phi, -1/phi, scale),\\\n",
    "                             vertex(1/phi, 0, phi, scale),\\\n",
    "                             vertex(-1/phi, 0, phi, scale),\\\n",
    "                             vertex(1/phi, 0, -phi, scale),\\\n",
    "                             vertex(-1/phi, 0, -phi, scale),\\\n",
    "                             vertex(phi, 1/phi, 0, scale),\\\n",
    "                             vertex(phi, -1/phi, 0, scale),\\\n",
    "                             vertex(-phi, 1/phi, 0, scale),\\\n",
    "                             vertex(-phi, -1/phi, 0, scale)]\n",
    "            edges = [(8,0), (8,1), (9,2), (9,4), (10,3), (10,5), (10,8), (11,6), (11,7),\\\n",
    "                     (11,9), (12,0), (12,2), (13,1), (13,4), (13,12), (14,3), (14,6), (15,5),\\\n",
    "                     (15,7), (15,14), (16,0), (16,3), (17,2), (17,6), (17,16), (18,1), (18,5),\\\n",
    "                     (19,4), (19,7), (19,18)]\n",
    "            faces = [(0,8,10,3,16),(0,12,13,1,8),(8,1,18,5,10),(10,5,15,14,3),(3,14,6,17,16),(16,17,2,12,0),\\\n",
    "                     (4,9,11,7,19),(4,13,12,2,9),(9,2,17,6,11),(11,6,14,15,7),(7,15,5,18,19),(19,18,1,13,4)]\n",
    "            edges = np.array(edges)\n",
    "            CoordVertices = np.array(CoordVertices)\n",
    "            faces = np.array(faces)\n",
    "        return CoordVertices, edges, faces\n",
    "\n",
    "    def coords(self):\n",
    "        # central atom = \"1st shell\"\n",
    "        c = [[0., 0., 0.]]\n",
    "        self.nAtoms = 1\n",
    "        self.nAtomsPerShell = [0]\n",
    "        self.nAtomsPerShell[0] = 1\n",
    "        indexVertexAtoms = []\n",
    "        indexEdgeAtoms = []\n",
    "        indexFaceAtoms = []\n",
    "        for i in range(1,self.nShell+1):\n",
    "            # vertices\n",
    "            nAtoms0 = self.nAtoms\n",
    "            cshell, E, F = self.MakeVertices(i)\n",
    "            self.nAtoms += self.nVertices\n",
    "            self.nAtomsPerShell.append(self.nVertices)\n",
    "            c.extend(cshell.tolist())\n",
    "            indexVertexAtoms.extend(range(nAtoms0,self.nAtoms))\n",
    "\n",
    "            # intermediate atoms on edges e\n",
    "            nAtoms0 = self.nAtoms\n",
    "            Rvv = RAB(cshell,E[0,0],E[0,1]) #distance between two vertex atoms\n",
    "            nAtomsOnEdges = int((Rvv+1e-6) / self.Rnn)-1\n",
    "            nIntervals = nAtomsOnEdges + 1\n",
    "            #print(\"nAtomsOnEdges = \",nAtomsOnEdges)\n",
    "            coordEdgeAt = []\n",
    "            for n in range(nAtomsOnEdges):\n",
    "                for e in E:\n",
    "                    a = e[0]\n",
    "                    b = e[1]\n",
    "                    coordEdgeAt.append(cshell[a]+vector(cshell,a,b)*(n+1) / nIntervals)\n",
    "            self.nAtomsPerShell[i] += nAtomsOnEdges * len(E) # number of edges x nAtomsOnEdges\n",
    "            self.nAtoms += nAtomsOnEdges * len(E)\n",
    "            c.extend(coordEdgeAt)\n",
    "            indexEdgeAtoms.extend(range(nAtoms0,self.nAtoms))\n",
    "            #print(c)\n",
    "            \n",
    "            # center of each pentagonal facet\n",
    "            nAtomsOnFaces = 0\n",
    "            nAtoms0 = self.nAtoms\n",
    "            coordFaceAt = []\n",
    "            for f in F:\n",
    "                nAtomsOnFaces += 1\n",
    "                coordCenterFace = centerOfGravity(cshell,f)\n",
    "                #print(\"coordCenterFace\",coordCenterFace)\n",
    "                self.nAtomsPerShell[i] += 1\n",
    "                coordFaceAt.append(coordCenterFace)\n",
    "                # atoms from the center of each pentagonal facet to each of its apex \n",
    "                nAtomsOnInternalRadius = i-1\n",
    "                nIntervals = nAtomsOnInternalRadius+1\n",
    "                print(f)\n",
    "                for indexApex,apex in enumerate(f):\n",
    "                    if (indexApex == len(f)-1):\n",
    "                        indexApexPlus1 = 0\n",
    "                    else:\n",
    "                        indexApexPlus1 = indexApex+1\n",
    "                    apexPlus1 = f[indexApexPlus1]\n",
    "                    for n in range(nAtomsOnInternalRadius):\n",
    "                        nAtomsOnFaces += 1\n",
    "                        coordAtomOnApex = coordCenterFace+vectorBetween2Points(coordCenterFace,cshell[apex])*(n+1) / nIntervals\n",
    "                        coordAtomOnApexPlus1 = coordCenterFace+vectorBetween2Points(coordCenterFace,cshell[apexPlus1])*(n+1) / nIntervals\n",
    "                        coordFaceAt.append(coordAtomOnApex)\n",
    "                        RbetweenRadialAtoms = Rbetween2Points(coordAtomOnApex,coordAtomOnApexPlus1)\n",
    "                        nAtomsBetweenRadialAtoms = int((RbetweenRadialAtoms+1e-6) / self.Rnn)-1\n",
    "                        nIntervalsRadial = nAtomsBetweenRadialAtoms + 1\n",
    "                        for k in range(nAtomsBetweenRadialAtoms):\n",
    "                            nAtomsOnFaces += 1\n",
    "                            coordFaceAt.append(coordAtomOnApex+vectorBetween2Points(coordAtomOnApex,coordAtomOnApexPlus1)*(k+1) / nIntervalsRadial)\n",
    "            self.nAtoms += nAtomsOnFaces\n",
    "            c.extend(coordFaceAt)\n",
    "            indexFaceAtoms.extend(range(nAtoms0,self.nAtoms))\n",
    "                \n",
    "            # now, facet atoms\n",
    "            # coordFaceAt = []\n",
    "            # nAtomsOnFaces = 0\n",
    "            # nAtoms0 = self.nAtoms\n",
    "            # for f in F:\n",
    "            #     print(\"f = \",f)\n",
    "            #     #nAtomsOnFaces,coordFaceAt = MakeFaceCoord(self.Rnn,f,cshell,nAtomsOnFaces,coordFaceAt)\n",
    "            # self.nAtomsPerShell[i] += nAtomsOnFaces\n",
    "            # self.nAtoms += nAtomsOnFaces\n",
    "            # c.extend(coordFaceAt)\n",
    "            # indexFaceAtoms.extend(range(nAtoms0,self.nAtoms))\n",
    "            \n",
    "        # print(indexVertexAtoms)\n",
    "        # print(indexEdgeAtoms)\n",
    "        # print(indexFaceAtoms)\n",
    "        return c,[indexVertexAtoms,indexEdgeAtoms,indexFaceAtoms]\n",
    "    \n",
    "    def prop(self):\n",
    "        print(self)\n",
    "        print(\"number of vertices = \",self.nVertices)\n",
    "        print(\"number of edges = \",self.nEdges)\n",
    "        print(\"number of faces = \",self.nFaces)\n",
    "        print(\"phi = \",self.phi)\n",
    "        print(f\"intershell factor = {self.interShellF:.2f}\")\n",
    "        print(f\"nearest neighbour distance = {self.Rnn:.2f} Å\")\n",
    "        print(f\"intershell distance = {self.interShellDistance:.2f} Å\")\n",
    "        print(f\"edge length = {self.edgeLength()*0.1:.2f} nm\")\n",
    "        print(f\"radius after volume = {RadiusSphereAfterV(self.volume()*1e-3):.2f} nm\")\n",
    "        print(f\"radius of the circumscribed sphere = {self.radiusCircumscribedSphere()*0.1:.2f} nm\")\n",
    "        print(f\"radius of the inscribed sphere = {self.radiusInscribedSphere()*0.1:.2f} nm\")\n",
    "        print(f\"area = {self.area()*1e-2:.1f} nm2\")\n",
    "        print(f\"volume = {self.volume()*1e-3:.1f} nm3\")\n",
    "        print(\"number of atoms per shell = \",self.nAtomsPerShellAnalytic())\n",
    "        print(\"total number of atoms = \",self.nAtomsAnalytic())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "68a17e2d-7d8a-466b-a51c-9e83d8446749",
   "metadata": {},
   "outputs": [],
   "source": [
    "# DD1 = regdodeca(2.7,1)\n",
    "# c1,i1 = DD1.coords()\n",
    "# n = 0\n",
    "# for i in range(1,DD1.nAtoms):#do not take into account first atom\n",
    "#     for j in range(1,i):\n",
    "#         if RAB(c1,i,j) < 2.8:\n",
    "#             n += 1\n",
    "#             print(f\"{n}  ({i-1},{j-1})\") #shift => do not take into account inner atom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "ac364beb-6ad3-4ce3-84b8-7c27a30edb12",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Regular dodecahedron with 5 shell(s) and Rnn = 2.7\n",
      "number of vertices =  20\n",
      "number of edges =  30\n",
      "number of faces =  12\n",
      "phi =  1.618033988749895\n",
      "intershell factor = 0.71\n",
      "nearest neighbour distance = 2.70 Å\n",
      "intershell distance = 3.78 Å\n",
      "edge length = 1.35 nm\n",
      "radius after volume = 0.69 nm\n",
      "radius of the circumscribed sphere = 1.89 nm\n",
      "radius of the inscribed sphere = 1.50 nm\n",
      "area = 37.6 nm2\n",
      "volume = 1.4 nm3\n",
      "number of atoms per shell =  [1, 32, 122, 272, 482, 752]\n",
      "total number of atoms =  1661\n",
      "[ 0  8 10  3 16]\n",
      "[ 0 12 13  1  8]\n",
      "[ 8  1 18  5 10]\n",
      "[10  5 15 14  3]\n",
      "[ 3 14  6 17 16]\n",
      "[16 17  2 12  0]\n",
      "[ 4  9 11  7 19]\n",
      "[ 4 13 12  2  9]\n",
      "[ 9  2 17  6 11]\n",
      "[11  6 14 15  7]\n",
      "[ 7 15  5 18 19]\n",
      "[19 18  1 13  4]\n",
      "[ 0  8 10  3 16]\n",
      "[ 0 12 13  1  8]\n",
      "[ 8  1 18  5 10]\n",
      "[10  5 15 14  3]\n",
      "[ 3 14  6 17 16]\n",
      "[16 17  2 12  0]\n",
      "[ 4  9 11  7 19]\n",
      "[ 4 13 12  2  9]\n",
      "[ 9  2 17  6 11]\n",
      "[11  6 14 15  7]\n",
      "[ 7 15  5 18 19]\n",
      "[19 18  1 13  4]\n",
      "[ 0  8 10  3 16]\n",
      "[ 0 12 13  1  8]\n",
      "[ 8  1 18  5 10]\n",
      "[10  5 15 14  3]\n",
      "[ 3 14  6 17 16]\n",
      "[16 17  2 12  0]\n",
      "[ 4  9 11  7 19]\n",
      "[ 4 13 12  2  9]\n",
      "[ 9  2 17  6 11]\n",
      "[11  6 14 15  7]\n",
      "[ 7 15  5 18 19]\n",
      "[19 18  1 13  4]\n",
      "[ 0  8 10  3 16]\n",
      "[ 0 12 13  1  8]\n",
      "[ 8  1 18  5 10]\n",
      "[10  5 15 14  3]\n",
      "[ 3 14  6 17 16]\n",
      "[16 17  2 12  0]\n",
      "[ 4  9 11  7 19]\n",
      "[ 4 13 12  2  9]\n",
      "[ 9  2 17  6 11]\n",
      "[11  6 14 15  7]\n",
      "[ 7 15  5 18 19]\n",
      "[19 18  1 13  4]\n",
      "[ 0  8 10  3 16]\n",
      "[ 0 12 13  1  8]\n",
      "[ 8  1 18  5 10]\n",
      "[10  5 15 14  3]\n",
      "[ 3 14  6 17 16]\n",
      "[16 17  2 12  0]\n",
      "[ 4  9 11  7 19]\n",
      "[ 4 13 12  2  9]\n",
      "[ 9  2 17  6 11]\n",
      "[11  6 14 15  7]\n",
      "[ 7 15  5 18 19]\n",
      "[19 18  1 13  4]\n",
      "1661\n",
      "[1, 32, 62, 92, 122, 152]\n",
      "Moments of inertia = 43141644.62 43141644.62 43141644.62 amu.Å2\n",
      "Moments of inertia / M = 708.84 708.84 708.84 amu.Å2\n",
      "Size of the ellipsoid = 11.91 11.91 11.91 nm\n"
     ]
    }
   ],
   "source": [
    "DD1 = regDD(2.7,5)\n",
    "DD1.prop()\n",
    "c1,i1 = DD1.coords()\n",
    "print(DD1.nAtoms)\n",
    "print(DD1.nAtomsPerShell)\n",
    "dd1 = Atoms('Au'*DD1.nAtoms, positions=c1)\n",
    "write(\"3DNP-coords/dd1.xyz\", dd1)\n",
    "dd1.moi = dd1.get_moments_of_inertia() # in amu*angstrom**2\n",
    "print(f\"Moments of inertia = {dd1.moi[0]:.2f} {dd1.moi[1]:.2f} {dd1.moi[2]:.2f} amu.Å2\")\n",
    "dd1.masses = dd1.get_masses()\n",
    "dd1.M = dd1.masses.sum()\n",
    "dd1.moiM = dd1.moi/ico1.M\n",
    "print(f\"Moments of inertia / M = {dd1.moiM[0]:.2f} {dd1.moiM[1]:.2f} {dd1.moiM[2]:.2f} amu.Å2\")\n",
    "dd1.dim = 2*np.sqrt(5*dd1.moiM)\n",
    "print(f\"Size of the ellipsoid = {dd1.dim[0]*0.1:.2f} {dd1.dim[1]*0.1:.2f} {dd1.dim[2]*0.1:.2f} nm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2735dac-8de7-4a2c-9eff-73ac68b7eace",
   "metadata": {},
   "source": [
    "# Geometry optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "affe58ad-044f-4d5e-beab-e343a77e19ec",
   "metadata": {},
   "source": [
    "# End"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "d18987e2-9935-4018-9011-a7f1a1bd2d1b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "**End at:** Friday 06 October 2023, 15:59:34  \n",
       "**Duration:** 00:00:38 359ms"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "<p style=\"text-align: center\"><img width=\"800px\" src=\"./config/svg/logoEnd.svg\" style=\"margin-left:auto; margin-right:auto\"/></p>"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "vID.end(cwd0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7c790eb-82a1-455f-8319-47999e7d0396",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 5
}