particles-simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Particles simulation\n",
    "Bibliotek Atooms dostarcza interfejs wysokiego poziomu do tworzenia symulacji opartych na cząsteczkach.\n",
    "W głównej mierze skupia się na dynamice molekularnej oraz symulacjach Monte Carlo, ale bibliotekę Atooms można sotosować także do innych celów."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import atooms\n",
    "import numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Podstawy\n",
    "Podstawowowym obiektem wykorzystywanym przez Atooms to ```Particle```. Obiekt ten posiada swoje włąsciwości, w szczególności\n",
    "położenie (jako tablica numpy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from atooms.system.particle import Particle\n",
    "particle = Particle(position=[2.0, 2.0])\n",
    "print(particle.position, type(particle.position))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cząstecteczki można definiować w wielu wymiarach, ale najbardziej skupimy się na dwóch wymiarach."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "particle = Particle(position=[3.0, 3.0, 2.0, 1.0, 0.0])\n",
    "print(len(particle.position))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Możemy również definiować takie właściwości jak prędkość, chemiczny perwiastek, masę czy promień."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "particle = Particle(position=[1.0, 1.0, 1.0], velocity=[1.0, 0.0, 0.0])\n",
    "particle.position += numpy.array([0.0, 1.0, 1.0])\n",
    "particle.velocity *= 2\n",
    "particle.mass = 10.\n",
    "particle.radius = 4. \n",
    "print(particle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Można także dodawać własne włąściwości np takie jak ładunek elektryczny, który będzie przydatny w fizyce."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "particle.charge = -1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Możemy też definiować komórki ```Cell``` które tworzą przestrzeń dla naszych cząsteczek"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from atooms.system.cell import Cell\n",
    "cell = Cell(side=[3.0, 3.0])\n",
    "print(cell.side, cell.volume)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Możemy też ustawić cząsteczki w środku naszej komórki, wykorzystując funckę ```fold```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell = Cell(side=[3.0, 3.0])\n",
    "particle = Particle(position=[4.0, 4.0])\n",
    "print(particle)\n",
    "particle.fold(cell)\n",
    "print(particle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## System\n",
    "Obiekt ```System``` to swoistego rodzaju środowisko, które zawiera wszystkie ważne inforamcje, fizyczne obiekty \n",
    "naszej symulacji. Można także definiować w nim rezerwuary takie jak ```thermostats```, ```barostats ```, które to \n",
    "obiekty przetrzymuj informację o temeraturze lub ciśnieniu. Obiekt ```System``` powinien też zawierać wszystkie \n",
    "inne klasy opisujące oddziaływania pomiedzy cząsteczkami.\n",
    "\n",
    "Poniżej znajduje się definicja prostego systemu z paroma cząsteczkami w komórce w 2D. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from atooms.system import System\n",
    "system = System(particle=[Particle() for i in range(10)],\n",
    "                cell=Cell([10.0, 10.0]))\n",
    "\n",
    "system.density = 1.2  # równoważne system.set_density(1.2)\n",
    "system.temperature = 1.7  # równoważne system.set_temperature(1.2)\n",
    "print(system.density, system.temperature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Według dokumentacji, temeratura ustawiona powyżje to temeratura kinetyczna i nie musi być równa temeraturze \n",
    "w obiekcie ```Thermostat```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from atooms.system import Thermostat\n",
    "system.thermostat = Thermostat(temperature=1.0)\n",
    "system.temperature = 1.5  # equivalent to system.set_temperature(1.2)\n",
    "print(system.temperature, system.thermostat.temperature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pomocnicze funcje"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_frame(particles, i):\n",
    "    global trajectory\n",
    "    if i == 0:\n",
    "        trajectory = ''\n",
    "    trajectory += str(len(particles)) + '\\ntitle\\n'\n",
    "    for particle in particles:\n",
    "        trajectory += ' '.join(['Ar',str(particle.position[0]),str(particle.position[1]),'0.0\\n'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from atooms.system import System\n",
    "global trajectory\n",
    "\n",
    "class SimpleBackend(object):\n",
    "\n",
    "    def __init__(self):\n",
    "        # There are no particles at the beginning\n",
    "        system = System()\n",
    "        assert len(system.particle) == 0\n",
    "\n",
    "        # Add particles\n",
    "        from atooms.system.particle import Particle\n",
    "        from random import random\n",
    "        L = 3\n",
    "        for i in range(10):\n",
    "            p = Particle(position=[L * random(), L * random()], velocity = [random(), random()], radius = 0.1)\n",
    "            system.particle.append(p)\n",
    "        self.system = system\n",
    "\n",
    "    def run(self, steps):\n",
    "        for i in range(steps):\n",
    "            for particle in self.system.particle:\n",
    "                particle.position += particle.velocity * 0.01\n",
    "            add_frame(self.system.particle, i)\n",
    "\n",
    "# The backend is created and wrapped by a simulation object.\n",
    "# Here we first call the run() method then run_until()\n",
    "from atooms.simulation import Simulation\n",
    "backend = SimpleBackend()\n",
    "simulation = Simulation(backend)\n",
    "simulation.run_until(300)\n",
    "assert simulation.current_step == 300"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import py3Dmol\n",
    "view = py3Dmol.view()\n",
    "view.addModelsAsFrames(trajectory,'xyz')\n",
    "view.animate({'loop': 'forward', 'reps': 1})\n",
    "view.setStyle({'sphere':{'radius': 0.1}})\n",
    "view.zoomTo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zadanie 1.\n",
    "\n",
    "Zaimplementować prostą funkcję detekcji kolizji.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PyMunk\n",
    "Przykłąd przedstawia zachowanie energi w układzie \"kulek\" Newtona"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.insert(1,'/opt/conda/lib/python3.7/site-packages/pymunk/libchipmunk.so')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import animation\n",
    "from IPython.display import HTML"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymunk\n",
    "from pymunk.vec2d import Vec2d\n",
    "import pymunk.matplotlib_util"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup_space():\n",
    "    space = pymunk.Space()\n",
    "    space.gravity = 0,-9820\n",
    "    space.damping = 0.99\n",
    "    return space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup_balls(space):\n",
    "    width = 600\n",
    "    height = 600\n",
    "    for x in range(-100,150,50):\n",
    "        x += width / 2\n",
    "        offset_y = height/2\n",
    "        mass = 10\n",
    "        radius = 25\n",
    "        moment = pymunk.moment_for_circle(mass, 0, radius, (0,0))\n",
    "        body = pymunk.Body(mass, moment)\n",
    "        body.position = x, -125+offset_y\n",
    "        body.start_position = Vec2d(body.position)\n",
    "        shape = pymunk.Circle(body, radius)\n",
    "        shape.elasticity = 0.9999999\n",
    "        space.add(body, shape)\n",
    "        pj = pymunk.PinJoint(space.static_body, body, (x, 125+offset_y), (0,0))\n",
    "        space.add(pj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.axes(xlim=(0, 600), ylim=(0, 600))\n",
    "ax.set_aspect(\"equal\")\n",
    "\n",
    "space = setup_space()\n",
    "setup_balls(space)\n",
    "\n",
    "space.shapes[1].body.apply_impulse_at_local_point((-12000,0))\n",
    "    \n",
    "def init():\n",
    "    space.debug_draw(pymunk.matplotlib_util.DrawOptions(ax))\n",
    "    return []\n",
    "\n",
    "def animate(dt):\n",
    "    #we run the animation with half speed intentionally to make it a little nicer to look at\n",
    "    for x in range(10):\n",
    "        space.step(1/50/10/2)\n",
    "    ax.clear()\n",
    "    ax.set_xlim((0, 600))\n",
    "    ax.set_ylim((0, 600))\n",
    "    ax.set_aspect(\"equal\")\n",
    "    space.debug_draw(pymunk.matplotlib_util.DrawOptions(ax))\n",
    "    return []\n",
    "\n",
    "frames = 105\n",
    "anim = animation.FuncAnimation(fig, animate, init_func=init,\n",
    "                               frames=frames, interval=20, blit=True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generowanie filmiku z symulacją\n",
    "HTML(anim.to_html5_video())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "\n",
    "def add_ball(space):\n",
    "    mass = 1\n",
    "    radius = 1\n",
    "    moment = pymunk.moment_for_circle(mass, 0, radius) # 1\n",
    "    body = pymunk.Body(mass, moment) # 2\n",
    "    x = random.randint(2, 18)\n",
    "    body.position = x, 15 # 3\n",
    "    shape = pymunk.Circle(body, radius) # 4\n",
    "    shape.elasticity = 0.95\n",
    "    shape.friction = 0.95\n",
    "    space.add(body, shape) # 5\n",
    "    return shape\n",
    "\n",
    "\n",
    "def add_static_square(space):\n",
    "    body = pymunk.Body(body_type = pymunk.Body.STATIC) # 1\n",
    "\n",
    "    l1 = pymunk.Segment(body, (1, 1), (19, 1), 0.0) # 2\n",
    "    l2 = pymunk.Segment(body, (1, 1), (1, 19), 0.0)\n",
    "    l3 = pymunk.Segment(body, (1, 19), (19, 19), 0.0)\n",
    "    l4 = pymunk.Segment(body, (19, 1), (19, 19), 0.0)\n",
    "    \n",
    "    l1.elasticity = 0.65\n",
    "    l2.elasticity = 0.65\n",
    "    l3.elasticity = 0.65\n",
    "    l4.elasticity = 0.65\n",
    "    \n",
    "    l1.friction = 0.99\n",
    "\n",
    "    space.add(l1, l2, l3, l4) # 3\n",
    "    return l1,l2, l3, l4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# \n",
    "fig = plt.figure()\n",
    "ax = plt.axes(xlim=(0, 20), ylim=(0, 20))\n",
    "ax.set_aspect(\"equal\")\n",
    "\n",
    "space = setup_space()\n",
    "screen = pymunk.matplotlib_util.DrawOptions(ax)\n",
    "\n",
    "add_ball(space)\n",
    "add_ball(space)\n",
    "add_ball(space)\n",
    "add_ball(space)\n",
    "add_static_square(space)\n",
    "    \n",
    "def init():\n",
    "    print(screen)\n",
    "    space.debug_draw(screen)\n",
    "    return []\n",
    "\n",
    "def animate(dt):\n",
    "    #we run the animation with half speed intentionally to make it a little nicer to look at\n",
    "    for x in range(10):\n",
    "        space.step(1/50/10/10)\n",
    "    ax.clear()\n",
    "    ax.set_xlim((0, 20))\n",
    "    ax.set_ylim((0, 20))\n",
    "    ax.set_aspect(\"equal\")\n",
    "    space.debug_draw(screen)\n",
    "    return []\n",
    "\n",
    "frames = 105\n",
    "anim = animation.FuncAnimation(fig, animate, init_func=init,\n",
    "                               frames=frames, interval=20, blit=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "HTML(anim.to_html5_video())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Zadanie 2.\n",
    "Wykoryżystując kwadrat powyżej, dodać na spodzie wasrstwę kulek o masie 1 (zostawiając miejsce po bokach), a następnie dodać drugą warstwę kulek o masie 3, tak aby się warstwy się stykały."
   ]
  }
 ],
 "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.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}