AtomicSolid.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div id=\"glowscript\" class=\"glowscript\"></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from vpython import *\n",
    "scene = canvas() # This is needed in Jupyter notebook and lab to make programs easily rerunnable\n",
    "\n",
    "# Bruce Sherwood\n",
    "\n",
    "N = 3 # N by N by N array of atoms\n",
    "# Surrounding the N**3 atoms is another layer of invisible fixed-position atoms\n",
    "# that provide stability to the lattice.\n",
    "k = 1\n",
    "m = 1\n",
    "spacing = 1\n",
    "atom_radius = 0.3*spacing\n",
    "L0 = spacing-1.8*atom_radius\n",
    "V0 = pi*(0.5*atom_radius)**2*L0 # initial volume of spring\n",
    "scene.center = 0.5*(N-1)*vector(1,1,1)\n",
    "dt = 0.04*(2*pi*sqrt(m/k))\n",
    "axes = [vector(1,0,0), vector(0,1,0), vector(0,0,1)]\n",
    "\n",
    "scene.caption= \"\"\"A model of a solid represented as atoms connected by interatomic bonds.\n",
    "\n",
    "Right button drag or Ctrl-drag to rotate \"camera\" to view scene.\n",
    "To zoom, drag with middle button or Alt/Option depressed, or use scroll wheel.\n",
    "  On a two-button mouse, middle is left + right.\n",
    "Touch screen: pinch/extend to zoom, swipe or two-finger rotate.\"\"\"\n",
    "\n",
    "class crystal:\n",
    "        \n",
    "    def __init__(self,  N, atom_radius, spacing, momentumRange ):\n",
    "        self.atoms = []\n",
    "        self.springs = []\n",
    "        \n",
    "        # Create (N+2)^3 atoms in a grid; the outermost atoms are fixed and invisible\n",
    "        for z in range(-1,N+1,1):\n",
    "            for y in range(-1,N+1,1):\n",
    "                for x in range(-1,N+1,1):\n",
    "                    atom = sphere()\n",
    "                    atom.pos = vector(x,y,z)*spacing\n",
    "                    atom.radius = atom_radius\n",
    "                    atom.color = vector(0,0.58,0.69)\n",
    "                    if 0 <= x < N and 0 <= y < N and 0 <= z < N:\n",
    "                        p = vec.random()\n",
    "                        atom.momentum = momentumRange*p\n",
    "                    else:\n",
    "                        atom.visible = False\n",
    "                        atom.momentum = vec(0,0,0)\n",
    "                    atom.index = len(self.atoms)\n",
    "                    self.atoms.append( atom )\n",
    "        for atom in self.atoms:\n",
    "            if atom.visible:\n",
    "                if atom.pos.x == 0:\n",
    "                    self.make_spring(self.atoms[atom.index-1], atom, False)\n",
    "                    self.make_spring(atom, self.atoms[atom.index+1], True)\n",
    "                elif atom.pos.x == N-1:\n",
    "                    self.make_spring(atom, self.atoms[atom.index+1], False)\n",
    "                else:\n",
    "                    self.make_spring(atom, self.atoms[atom.index+1], True)\n",
    "\n",
    "                if atom.pos.y == 0:\n",
    "                    self.make_spring(self.atoms[atom.index-(N+2)], atom, False)\n",
    "                    self.make_spring(atom, self.atoms[atom.index+(N+2)], True)\n",
    "                elif atom.pos.y == N-1:\n",
    "                    self.make_spring(atom, self.atoms[atom.index+(N+2)], False)\n",
    "                else:\n",
    "                    self.make_spring(atom, self.atoms[atom.index+(N+2)], True)\n",
    "                    \n",
    "                if atom.pos.z == 0:\n",
    "                    self.make_spring(self.atoms[atom.index-(N+2)**2], atom, False)\n",
    "                    self.make_spring(atom, self.atoms[atom.index+(N+2)**2], True)\n",
    "                elif atom.pos.z == N-1:\n",
    "                    self.make_spring(atom, self.atoms[atom.index+(N+2)**2], False)\n",
    "                else:\n",
    "                    self.make_spring(atom, self.atoms[atom.index+(N+2)**2], True)\n",
    "    \n",
    "    # Create a grid of springs linking each atom to the adjacent atoms\n",
    "    # in each dimension, or to invisible motionless atoms\n",
    "    def make_spring(self, start, end, visible):\n",
    "        spring = helix()\n",
    "        spring.pos = start.pos\n",
    "        spring.axis = end.pos-start.pos\n",
    "        spring.visible = visible\n",
    "        spring.thickness = 0.05\n",
    "        spring.radius = 0.5*atom_radius\n",
    "        spring.length = spacing\n",
    "        spring.start = start\n",
    "        spring.end = end\n",
    "        spring.color = color.orange\n",
    "        self.springs.append(spring)\n",
    "\n",
    "c = crystal(N, atom_radius, spacing, 0.1*spacing*sqrt(k/m))\n",
    "\n",
    "while True:\n",
    "    rate(60)\n",
    "    for atom in c.atoms:\n",
    "        if atom.visible:\n",
    "            atom.pos = atom.pos + atom.momentum/m*dt\n",
    "    for spring in c.springs:\n",
    "        spring.axis = spring.end.pos - spring.start.pos\n",
    "        L = mag(spring.axis)\n",
    "        spring.axis = spring.axis.norm()\n",
    "        spring.pos = spring.start.pos+0.5*atom_radius*spring.axis\n",
    "        Ls = L-atom_radius\n",
    "        spring.length = Ls\n",
    "        Fdt = spring.axis * (k*dt * (1-spacing/L))\n",
    "        if spring.start.visible:\n",
    "            spring.start.momentum = spring.start.momentum + Fdt\n",
    "        if spring.end.visible:\n",
    "            spring.end.momentum = spring.end.momentum - Fdt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}