pkdtree test

pkdtree.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from MDAnalysis.lib import distances"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generate test particles "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = 5\n",
    "box = np.array([[d, .5 *d , .5*d],\n",
    "                [0, .5 * np.sqrt(3) *d, 1/6*np.sqrt(3)*d],\n",
    "                [0, 0, 1/3 * np.sqrt(6) * d]], dtype=np.float32).T\n",
    "#box = np.eye(3).astype(np.float32) * d\n",
    "\n",
    "a, b, c = box\n",
    "\n",
    "num = 500\n",
    "pos = np.random.uniform(0, 1, size=(num, 3)).astype(np.float32)\n",
    "pos = pos[:, 0][:, np.newaxis] * a + pos[:, 1][:, np.newaxis] * b "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualize Test particles "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = []\n",
    "for i in range(-1, 3):\n",
    "    for j in range(-1, 3):\n",
    "        grid.append(i * a + b * j)\n",
    "grid = np.asarray(grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEYCAYAAACHjumMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt4VfWd7/H3N8nG7ARIUPBCAEFULhIIGhABqzV26EWR\n8TijPaO19ZwybcdWZ1pmtNbWsY9VD6dTevdhOr34FK1ToKlRWy8R7SHeCAYJqBFBjASsEZtdYO/A\nTvI7fyTZcgmQZK+VtdfO5/U8PJC1117r+8te+8Nv3X7LnHOIiPghJ+gCRCR7KWBExDcKGBHxjQJG\nRHyjgBER3yhgRMQ3ChgR8Y0CRkR8o4AREd/kBbHSkSNHuvHjxwexahHxwPr16993zo063nyBBMz4\n8eOpra0NYtUi4gEze7s382kXSUR8o4AREd8oYETEN4Ecg+lJMplkx44dtLa2Bl1KqOTn5zNmzBgi\nkUjQpYgcIWMCZseOHQwbNozx48djZkGXEwrOOXbv3s2OHTuYMGFC0OWIHCFjdpFaW1s56aSTFC59\nYGacdNJJ6vVJxsqYHgygcOmjeDzOnj17aGlpYdmyZVRUVFBaWhp0WX1SX19PdXU1sViMoqKi0LUh\n7PWDv23IqICR3ovH48RiMbqHPI3FYlRVVQGEZgOvr6+nqqqKZDIJhK8NYa8f/G9DxuwihcV9993H\n/ffff8T07du3M23atAGrY8+ePRw+nnIymaS6unrAakhXdXV1asPuFqY2hL1+8L8Noe3BBNE1bWtr\n4wtf+IKv6+it9vb2HqfHYrEBrqT/jlZrWNoQ9vrB/zaEMmD86tZ9+9vfZsWKFYwdO5aRI0dy3nnn\n8cgjjzB37lxqampYuHAhe/bsYejQoXzta19j/fr13HDDDRQUFDB//nxP2tZbubm5PYZMUVHRgNaR\njqKioh435LC0Iez1g/9tCOUukh/dutraWlatWkVdXR2rV68+5F6plpYWnn32Wb761a8e8p7Pfe5z\n/OAHP+D555/v93r7a9iwYUccFI9EIlRUVAx4Lf1VUVFxxPU7YWpD2OsH/9sQyh6MH926tWvXcsUV\nVxCNRgG4/PLLU69dffXVPa6rpaWFiy66CIDrrruOP/zhD/1ef18VFBQAncdigFCeweiuNaxnYcJe\nP/jfhlAGjB/dumM9gK6wsLDH+YM+rV5QUEBBQQEffPABN998c6C19FdpaWmovpCHC3v94G8bQrmL\n5Ee3bv78+VRVVdHa2srevXt59NFHjzl/cXExRUVFrF27FoAVK1b0e90i2SqUPRg/unWzZs1i4cKF\nzJgxg9NPP53y8vLj9oh+8YtfpA7yLliwoN/rFslWFsSzqcvLy93hA0699tprTJkyZcBrOdjevXsZ\nOnQo8Xicj3zkIyxfvpxzzz030Jp6IxN+dzK4mNl651z58ebzpAdjZv8M/G/AAfXA55xzobtBZvHi\nxbz66qu0trZy/fXXhyJcRDJZ2gFjZiXAV4CpzrmEmf03cA3wy3SXPdAeeOCBoEsQySpeHeTNA6Jm\nlgcUADs9Wq6IhFjaAeOcawL+L9AI7AJizrkn0l2uiIRf2gFjZiOAK4AJwGig0Myu7WG+xWZWa2a1\nzc3N6a5WRELAi12kS4G3nHPNzrkksBqYe/hMzrnlzrly51z5qFHHfZyKiGQBLwKmEZhjZgXWeWlr\nBfCaB8sNncrKSl599dXUz9/85jd56qmnjjr/M888w2WXXTYQpYkEIu2zSM65F81sJfAy0AbUAcvT\nXe7xVNY1sfTxBna2JBhdHGXJgkksmlni92qPqq2tjcrKSi677DKmTp0KwJ133hlYPSKZwJOzSM65\nbznnJjvnpjnnrnPO7fdiuUdTWdfEravraWpJ4ICmlgS3rq6nsq4preVu376dyZMnc/311zN9+nSu\nuuoq4vE4d955J7NmzWLatGksXrw4dd/SxRdfzNe//nUuuugi7r33Xh5++GGWLFlCWVkZW7du5bOf\n/SwrV64EYN26dcydO5cZM2Ywe/bs1E2K3fbt28cNN9zArFmzmDlzJr///e8B2Lx5M7Nnz6asrIzp\n06ezZcuWtNooMpBCeS/S0scbSCQPHQslkWxn6eMNaS+7oaGBxYsXs3HjRoYPH85PfvITbrzxRtat\nW8emTZtIJBI88sgjqfm7h3K47bbbWLhwIUuXLmXDhg1MnDgxNc+BAwe4+uqr+f73v88rr7zCU089\nlbpru9tdd93FJZdcwrp161izZg1Llixh37593Hfffdx0001s2LCB2tpaxowZk3YbRQZKKANmZ0ui\nT9P7YuzYscybNw+Aa6+9lrVr17JmzRrOP/98SktLefrpp9m8eXNq/p6GcjhcQ0MDp512GrNmzQJg\n+PDh5OUdunf6xBNPcM8991BWVsbFF19Ma2srjY2NXHDBBXznO9/h3nvv5e233z4imEQyWShvdhxd\nHKWphzAZXZz+l+/wIRjMjC996UvU1tYyduxY7rjjjkMeE9LTUA6H683QDs45Vq1axaRJkw6ZPmXK\nFM4//3weffRRFixYwM9+9jMuueSSPrRIJDih7MEsWTCJaCT3kGnRSC5LFkw6yjt6r7GxMTVC3YMP\nPpgaCnPkyJHs3bs3dUylJ8OGDTvi2ArA5MmT2blzJ+vWrQM6B4lqa2s7ZJ4FCxbwwx/+MHV8p66u\nDoBt27Zxxhln8JWvfIWFCxeycePGtNsoMlBCGTCLZpZw95WllBRHMaCkOMrdV5Z6chZpypQp/OpX\nv2L69Ol88MEHfPGLX+Tzn/88paWlLFq0KLWb05NrrrmGpUuXMnPmTLZu3ZqaPmTIEB566CG+/OUv\nM2PGDD72sY8d8bC022+/nWQyyfTp05k2bRq33347AA899BDTpk2jrKyM119/nc985jNpt1FkoGi4\nhoNs376dyy67jE2bNgVaR19lwu9OBpfeDtcQyh6MiISDAuYg48ePD13vRSSTZVTABLG7Fnb6nUkm\ny5iAyc/PZ/fu3frC9IFzjt27d5Ofnx90KSI9ypjrYMaMGcOOHTvQUA59k5+fr6t7JWNlTMBEIhEm\nTJgQdBki4qGM2UUSkeyjgBER3yhgRMQ3ChgR8Y0CRkR8o4AREd8oYETENwoYEfGNAkZEfKOAERHf\nKGBExDcKGBHxjQJGRHyjgBER33gSMGZWbGYrzex1M3vNzC7wYrkiEm5ejQfzfeCPzrmrzGwIUODR\nckUkxNIOGDMbDnwE+CyAc+4AcCDd5YpI+Hmxi3QG0Az8wszqzOxnZnbE81TNbLGZ1ZpZrYbFFBkc\nvAiYPOBc4KfOuZnAPuCWw2dyzi13zpU758pHjRrlwWpFJNN5ETA7gB3OuRe7fl5JZ+CIyCCXdsA4\n594F3jGz7ifPVwCvprtcEQk/r84ifRlY0XUGaRvwOY+WKyIh5knAOOc2AMd9ELaIDC66kldEfKOA\nERHfKGBExDcKGBHxjQJGRHyjgBER3yhgRMQ3ChgR8Y0CRkR8o4AREd8oYETENwoYEfGNAkZEfKOA\nERHfKGBExDcKGBHxjQJGRHyjgBER3yhgRMQ3ChgR8Y1XTxWQHtTX11NdXU0sFqOoqIiKigpKS0uD\nLqvX+lN/ZV0TSx9vYGdLgtHFUZYsmMSimSW9fr2/KuuauOPhzbQkkgCMKIjwrcvPYWLeB6H+DCDc\n25E55wZ8peXl5a62tnbA1zuQ6uvrqaqqIplMpqZFIhEuv/zyUGwcB9efSCTo6OggEomwYMECpk6d\n2uN7Htu4i28/+iqtyfZDphdHIyxZMBngiNfzI7lcPuM01m7ZzbuxBKcWRbnxo2fyyemn8djGXfxo\nzZup6fPPOom1W3azK5Yg14x251J/9yQvB+YM2cFY3icajZKTkxOqzwAydzsys/XOueM+SUQB45Nl\ny5YRi8WOmF5UVMTNN98cQEV9013/pk2bWLlyZdDlpK2srIxFixYB4fkMIHO3o94GjHaRfNLTRnGs\n6ZkmFouxte1E/lD3FrlDT2L4nKtSrw3JzWHhjNFMH1ucmrbxnRZWvrwjiFKPqaN1L7G1KygsLExN\nC8tnAOHfjhQwPikqKjrq/zxhsCtSwv9rKWTf9lcYfu7lDD/v8kNe31wcZfmNlwCdxz9+vLqe4efN\nDKLUY9pb/xQAwyfNSU0Ly2cA4d+OPDuLZGa5ZlZnZo94tcwwq6ioIBKJHDItEolQUVERUEV9s75t\nDHveXA/tbRRMmnvE600tCSbc8ijz7nmaf6/aTOKw4y6ZIt5QQ+7wUbw+8kK2tp0Yqs8Awr8dedmD\nuQl4DRju4TJDq/sAXBiP/lfWNdEcb+/8cg49iSGjJ/U4n6MzaDJVx/44ie11DJv5KTosj7r2sSy5\nfGIoPoNuYd6OwKOAMbMxwKeAu4B/8WKZ2aC0tDQ0G0K3yrombl1dT8eBBK1vvczQGQswC+flUomt\nL3X1wOYDsLcjErrPA8K5HXXzastZBvwr0HG0GcxssZnVmlltc3OzR6sVry19vIFEsp3E1lpc2wEK\nJs0LuqRjKimOMqIg0uNr+xpqyB16IieUdPbARhdHB7I0wYOAMbPLgPecc+uPNZ9zbrlzrtw5Vz5q\n1Kh0Vys+2dm1yxNvqCGnsJgTSqYEXNGx7dvfxl/iySOmdxxI0LptPQWT5mGWQzSSy5IFPe/qiX+8\n6MHMAxaa2XbgN8AlZvZrD5YrARhdHKUj2Upi2zoKzp6L5eQOyHpHFEQoiPR9c+y+cvdwqR7Y2XMp\nKY5y95WlnlwxLH2TdsA45251zo1xzo0HrgGeds5dm3ZlEoglCybR0bgBl9w/oLtHf4knGVF4gmfL\ni7/xHDkFxRRPKPXsdgTpu3AevRPfLJpZwsR9m8iJDid/7LQ+vz/XrF/rzbEPd8/S1ZFsJbF1HQVn\nX0Bre+dxJQmGpwHjnHvGOXeZl8uUgdXa2krd2mpGTpvfr92j4dH+nZjscJ2nvb3Q+tbLuGRrqgfm\nVXBJ36kHI4d48skn2bNnD0Mn92/3qKcDrgNtX0NNZw9sXOepXZ09Co4CRg5x23/8Jzn5Q2k75Zyg\nS+kX15Yk8eZLFJw1J9UD09mj4ChgJOW3L73FpuefJnrmHCw3nLepJbbX4Q4kUrtHIwoiOsAbIAWM\npHzrpw/h9u/r8d6jsIg31JBzQiH5p08H4FuXh7Mnli0UMJLyzstrsCFRouMz767o3nDtSRJbXiB6\n1hwsN0I0kqPeS8AUMAJAMpmkdcvzFJx5PpbX86X3ma717Y10HNQDu/vK6QFXJAoYAeDZZ5+lLbGH\nYVPmB11Kv8UbalI9sEgO6r1kAAWMALBq1SoKCwsZfuZ5QZfSL66jnfiWF4ieORvLG8LSvysLuiRB\nASNAe3s7q1evpmzuJewnpLtHjfV0JP5KYdfZI/VeMkM4z0WKp9auXct7771HtDi8xyziDTVY5ATy\nJ5zb79sVxHvqwQgrV65kSH4+HWPCuVvRuXv0PNEzZpETyT/qY0xk4ClgBrmOjg4eeOi3RMbNJGfI\n8S+pz8vJvN7B/qbX6NjXkrq4Tj2YzKGAGeReeOEFPmj+M9FeDs3Q3uGIRg69CTIayaU4Gtyxm3hD\nDZY3hOjEzsf0qAeTORQwg9zKlSshN4/oxNm9mt8Bd19ZSklxFIPUYE53LDzniOAZCM51EG+oIX/C\nuakeWIlubswYCphBzDnHqlWriE44l5wTCnr1HjOOeLY0fDiW70A7sLOB9r0fpHaPNDRmZtFZpEGs\ntraWxsZGxv3t3/f6PcaHjyppakmw5LevgEGyPZjdkn0NNZCbR8GZsynpCjydos4c6sEMYitXriQv\nL4+7bvpsr9/TcViOJDtcYOHinCPeUEN0/EwKhw5XuGQgBcwg5Zxj5cqVXHrppVx70TmBHqTtrwPv\nbqH9r80UTJpHItmuoTEzkAJmkNqwYQPbtm1j4uwK5t3z9FFH589k8YbnICeX6JnnAxoaMxPpGMwg\ntWrVKiwnl4f/Mpqc/b3/YkZyjOTh+0kB6N49yj99BrnRYYCGxsxE6sEMQs45fvHrBzlh7DRyCor6\n9N4OoDgaweh8EkBQku+9RVvLLp09ynAKmEFo8+bN7Hx7W7+ee9Te4djT2gZAkNez7WuoAcuh4Kw5\n5JrpwWoZSgEzyFTWNfGJm5cCRsHZF/RrGe3O4fDuMSN9ldo9GjeNoUUn8t2/n6FwyVAKmEGksq6J\nW1fX8+dXnuGEseeQWzgi6JL6pX13I20f7GB02UfVc8lwOsg7iCx9vIG/vrud5PuNjLj0H4Mup99y\ntr+EmVFz3y2ceuqpQZcjx5B2D8bMxprZGjN7zcw2m9lNXhQm3mtqSXSe2oV+7x5lgndfeYb58+cr\nXELAi12kNuCrzrkpwBzgn8xsqgfLFY/lmrGvoYYTSqaQN2xk0OX0S3L3DpLvv81VV10VdCnSC2kH\njHNul3Pu5a5/7wFeA7RTnIFaP9hJ8r1tFJwd4ucevdHZA7vyyisDrkR6w9ODvGY2HpgJvOjlcsUb\nkcaXAPp1ejpTxBtqGDZuKmPGjAm6FOkFzwLGzIYCq4CbnXN/7eH1xWZWa2a1zc3NXq1W+iDyzkvk\njz6bvKKTgy6lX5J/2cWBP28lb2J4jx8NNp4EjJlF6AyXFc651T3N45xb7pwrd86Vjxo1yovVSi9V\n1jVx3r89yJubX2HEORemrsQN28iS8TdqAJgw65KAK5HeSvs0tZkZ8F/Aa865/0i/JPFS6tqXddUA\n2IQ5xBLJzovkgr+lqE/iDc8x5NSz+MQFpUGXIr3kRQ9mHnAdcImZbej680kPlise6B5pLv7Gc0RO\nPoPIiNPClisAtMXe48CuNyiYNI9HN+4KuhzppbR7MM65tXQOdCYZaGdLgrY977O/6TWKL7wu6HL6\nrfvsUcGkufwlHr6hJQYr3SqQ5UYXR4m/8TyQmWePejvQVbyhhsjJE4iMGO1zReIlBUyWW7JgUueX\nc+Q4Iif1fGq3+zlCQXRD97d1HHee7h5Yd0CGcfS9wUoBk+Vadjez/53Nx+y9fPr8sWy/51N87+qB\nf7Jjb55E0N0DK5w0j0iOccfCc/wuSzyimx2z3DeW/RxwxwyYFS82suKFRkYXRymI5BBPHr9XMZDa\nt75AZOQ4xk88WwN7h4wCJotV1jWxe9OfyDtxDJGRpx91vu6Bo5oycExbF29hz/aN3H777fz7Lbr+\nJWy0i5TFvrP6JVob6ymYNBcL21V1XfY2PIdzTjc3hpQCJottrX0GXAeFGXj2qLfiDc8x7JRxTJs2\nLehSpB8UMFmqsq6JREMNecWnEjn5jKDL6Zf2eIzWxo3YGXNC2wMb7BQwWeruyloSb2+gYNK80H45\n41teBNcR6uElBjsFTJbaVvssdLT3+csZyc2cMIo31JBXdApDTpkYdCnSTwqYLPSNynr2vVFD7rBR\nDDnt7F6/r6Q4ytWzxqYuvAtSe+teWt9+JdUDq6xrCrok6QcFTJb5RmU99z/7Oom3Xu7T2aNlV5ex\nZMEkVq1voj3IBx51Sbz5InS0pa7f0XOnw0nXwWSZFS80ktj6ErS3UTBpfq/fd/NDG8gxyICnwgKd\nu0cH98D03OlwUg8myzg6n3qYO/RETijp26NUMyVcOvbHj+iB6bnT4aSAyTIdBxK0bltPwdlzMQvn\nx/thD+zD63f03OlwCucWKEflGutwbQcycmiG3oo3PNfVA5sMwLVzxun+o5DSMZgsc/pf69lVUMwJ\nY8LzaKocg+H5EWKJJMNy22jcVsvQ6R9jzIhC3dwYcgqYLBKPx1m/tpriqRdjOblBl9NrHQ5aEklK\niqPMzdlCfdsBHv7u17j44ouDLk3SpIDJInf+9AFaE3GGnxnOK1+bWhLc9/AKho84iQsvvDDocsQD\nOgaTRf7r178hJzqc/LHe3xhYUhylxOczOR3J/ex7cx05E2aTmxueHpgcnQImS+zfv5/drz5HwVlz\nsFzvO6YfnTyKJQsmEY3498VvfetlXLKVIWfO1ZW7WUIBkyXuvO83uAOJHu89KhySS0lxFKOzJ3Lt\nnHF9Xv6a15tZNLOE/3Gefwdc4w015OQPI39cqa7czRI6BpMFKuua+NHPf42dUEj++BlHvB7JzaHm\nsNHg1rze3KcR7Ha2JKisa2LVen96Fq4tSfzNlzrvPcrN05W7WUI9mCxw76Ob2NvwPAVnnY/lHjni\nfixx5HOEliyY1Kc7p4uiEf75oQ29GqS7PxJvb8AdiKcGx9KVu9lBAZMF3tr4Ih379x314rqevqyL\nZpZQOKR3HdhIjvHX1qSvT4SMv16T6oFFI7m6cjdLaBcpC7jtL2JDokTHz+zx9aaWBBNvfYx25ygp\njvLRyaNY83ozLT30bA5nBkmfb1Jy7UkSW56n4MzZ5ORFuPvKUl1clyU8CRgz+zjwfSAX+Jlz7h4v\nluun+vp6qquricViFBUVUVFRQWlpuB6qXl9fzxNPPEFs07NEJ87C8oYcdd7uIRiaWhL8+oXGXq9j\nIEZuaH17Y6oH9i8XnBiqcMmG7chPae8imVku8GPgE8BU4NNmltHXqdfX11NVVUUsFgMgFotRVVVF\nfX19wJX1XncbNm7cyP7EvnAP7P3Gc9iQKNMnjmHPpjWh+RyyYTvymxc9mNnAm865bQBm9hvgCuBV\nD5bti+rqapLJQ3cPkskk1dXVofnfp7sN999/PwB7ayvZs+GPAVfVP/t3vsbpZ53DhYXvkkwSms8h\nG7Yjv3kRMCXAOwf9vAM4//CZzGwxsBhg3Li+X4fhpe7/cXo7PRMdXmuhS7DnQIevB2K9loNjqB3g\n5FNPpmLOh4+tDcvnkA3bkd+8CJieznUesZ0755YDywHKy8sD/R4UFRX1uBEUFRUFUE3/dLfhjjvu\nSE3b2nYiLybHsd+nY/fRSC53X9n5P/PSxxv69SRIA753dRmLZpawbNmyUH8O2bAd+c2L09Q7gLEH\n/TwG2OnBcn1TUVFBJHLo9SKRSISKioqAKuq7ntowObqHlf9wBsuuLkvdN9Q9gHf3391X8hZHj7xe\nJpJrXDtnXOqq3xEFEYqjkdQVwN1ndxbNLKHmlktYdnUZkZzeX0tjwD8cNLZL2D+HsNc/EMyleZrA\nzPKAN4AKoAlYB/xP59zmo72nvLzc1dbWprXedGXD0f9021BZ18TSxxvY2ZJgdHG0X2OvHL6M7lPg\nO1sSFBdEcK7zQr+jLT/sn0PY6+8vM1vvnCs/7nzpBkzXyj4JLKPzNPXPnXN3HWv+TAgYEem/3gaM\nJzvrzrnHgMe8WJaIZA/dKiAivlHAiIhvFDAi4hsFjIj4RgEjIr5RwIiIbxQwIuIbBYyI+EYBIyK+\nUcCIiG8UMCLiGwWMiPhGASMivlHAiIhvFDAi4hsFjIj4RgEjIr5RwIiIbxQwIuIbBYyI+EYBIyK+\nUcCIiG8UMCLiGwWMiPhGASMivkkrYMxsqZm9bmYbzex3ZlbsVWEiEn7p9mCeBKY556YDbwC3pl+S\niGSLtALGOfeEc66t68cXgDHplyQi2cLLYzA3AH842otmttjMas2strm52cPVikimyjveDGb2FHBq\nDy/d5pz7fdc8twFtwIqjLcc5txxYDlBeXu76Va2IhMpxA8Y5d+mxXjez64HLgArnnIJDRFKOGzDH\nYmYfB/4NuMg5F/emJBHJFukeg/kRMAx40sw2mNl9HtQkIlkirR6Mc+5MrwoRkeyjK3lFxDcKGBHx\njQJGRHyjgBER3yhgRMQ3ChgR8Y0CRkR8o4AREd8oYETENwoYEfGNAkZEfKOAERHfKGBExDcKGBHx\njQJGRHyjgBER3yhgRMQ3ChgR8Y0CRkR8o4AREd8oYETENwoYEfGNAkZEfKOAERHfKGBExDcKGBHx\njScBY2ZfMzNnZiO9WJ6IZIe0A8bMxgIfAxrTL0dEsokXPZjvAf8KOA+WJSJZJK2AMbOFQJNz7pVe\nzLvYzGrNrLa5uTmd1YpISOQdbwYzewo4tYeXbgO+DvxNb1bknFsOLAcoLy9Xb0dkEDhuwDjnLu1p\nupmVAhOAV8wMYAzwspnNds69m25h9fX1VFdXE4vFKCoqoqKigtLS0nQXO6DC3oaw1w/Z0YYwO27A\nHI1zrh44uftnM9sOlDvn3k+3qPr6eqqqqkgmkwDEYjGqqqoAQrNxhL0NYa8fsqMNYZeR18FUV1en\nNopuyWSS6urqgCrqu7C3Iez1Q3a0Iez63YM5nHNuvFfLisVifZqeicLehrDXD9nRhrDLyB5MUVFR\nn6ZnorC3Iez1Q3a0IewyMmAqKiqIRCKHTItEIlRUVARUUd+FvQ1hrx+yow1h59kukpe6D8CF+eh/\n2NsQ9vohO9oQdubcwF+SUl5e7mprawd8vSLiDTNb75wrP958GbmLJCLZQQEjIr5RwIiIbxQwIuIb\nBYyI+EYBIyK+UcCIiG8CuQ7GzJqBtwd8xT0bCaR9B3jA1Ibghb1+6FsbTnfOjTreTIEETCYxs9re\nXDCUydSG4IW9fvCnDdpFEhHfKGBExDcKmK5xgkNObQhe2OsHH9ow6I/BiIh/1IMREd8oYETENwoY\nwMzuMLMmM9vQ9eeTQdfUW2b2cTNrMLM3zeyWoOvpKzPbbmb1Xb/3UAwSZGY/N7P3zGzTQdNONLMn\nzWxL198jgqzxeI7SBs+/BwqYD33POVfW9eexoIvpDTPLBX4MfAKYCnzazKYGW1W/fLTr9x6W60h+\nCXz8sGm3ANXOubOA6q6fM9kvObIN4PH3QAETbrOBN51z25xzB4DfAFcEXFPWc879CfjgsMlXAL/q\n+vevgEUDWlQfHaUNnlPAfOhGM9vY1XXM6O7tQUqAdw76eUfXtDBxwBNmtt7MFgddTBpOcc7tAuj6\n++TjzJ+pPP0eDJqAMbOnzGxTD3+uAH4KTATKgF3AdwMttvesh2lhu+5gnnPuXDp38/7JzD4SdEGD\nmOffg4x8qoAfjvaM7cOZ2X8Cj/hcjld2AGMP+nkMsDOgWvrFObez6+/3zOx3dO72/SnYqvrlz2Z2\nmnNul5mdBrwXdEF95Zz7c/e/vfoeDJoezLF0bRDd/hbYdLR5M8w64Cwzm2BmQ4BrgIcDrqnXzKzQ\nzIZ1/xuk6KttAAAAnUlEQVT4G8Lzuz/cw8D1Xf++Hvh9gLX0ix/fg0HTgzmO/2NmZXTuXmwH/jHY\ncnrHOddmZjcCjwO5wM+dc5sDLqsvTgF+Z2bQuS0+4Jz7Y7AlHZ+ZPQhcDIw0sx3At4B7gP82s/8F\nNAJ/F1yFx3eUNlzs9fdAtwqIiG+0iyQivlHAiIhvFDAi4hsFjIj4RgEjIr5RwIiIbxQwIuKb/w/m\nCkNRbw1f2gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6c9520d470>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=plt.figaspect(1))\n",
    "ax.scatter(grid[:, 0], grid[:, 1], color='#808080', label='grid')\n",
    "ax.scatter(pos[:, 0], pos[:, 1], label='particles')\n",
    "ax.plot([0, a[0], a[0] + b[0], b[0], 0], [0, 0, b[1], b[1], 0], color='black', zorder=4)\n",
    "ax.legend()\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Search for neighbors "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "pkd = mda.lib.pkdtree.PeriodicKDTree(box)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "pkd.set_coords(pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "center = np.asarray([d, d, 0], dtype=np.float32)\n",
    "center = (.5 * a + .5 * b) * 5\n",
    "\n",
    "pkd.search(center=center, radius=1)\n",
    "neighbors = pkd.get_indices()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAEYCAYAAACHjumMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8VPWd//HXJ5MJmQkkQYIXEhBE5WISEg2IgJcS2/Si\nyFq3tlutl9/K2taquxVX67ay9lEvpd2irZYH62r1IVpdoJFordaIdhNvBIIJqBGBGBOwRjQjMBMy\nyXx/fyQZISRkJnNOZs7k83w8eEDOnDnnczJn3nzPOd/zPWKMQSml7JAS7wKUUslLA0YpZRsNGKWU\nbTRglFK20YBRStlGA0YpZRsNGKWUbTRglFK20YBRStkmNR4rzcnJMZMnT47HqpVSFti0adMnxpjx\ng80Xl4CZPHkyNTU18Vi1UsoCIvJBJPPpIZJSyjYaMEop22jAKKVsE5dzMP0JBoM0NzfT3t4e71Ic\nJT09nby8PNxud7xLUeoICRMwzc3NjBkzhsmTJyMi8S7HEYwx7N27l+bmZqZMmRLvcpQ6QsIcIrW3\ntzNu3DgNlyiICOPGjdNWn0pYCdOCATRcouT3+9m3bx9tbW2sWLGC0tJSCgoK4l1WVOrr66msrMTn\n85GVleW4bXB6/WDvNiRUwKjI+f1+fD4fvUOe+nw+KioqAByzg9fX11NRUUEwGASctw1Orx/s34aE\nOURyipUrV/Loo48eMb2xsZH8/Pxhq2Pfvn30HU85GAxSWVk5bDXEqrKyMrxj93LSNji9frB/Gxzb\ngolH07Szs5Nrr73W1nVEqqurq9/pPp9vmCsZuoFqdco2OL1+sH8bHBkwdjXrfv7zn7N69WomTpxI\nTk4OZ5xxBs888wzz5s2jurqaRYsWsW/fPkaPHs1NN93Epk2buPrqq/F6vSxYsMCSbYuUy+XqN2Sy\nsrKGtY5YZGVl9bsjO2UbnF4/2L8NjjxEsqNZV1NTw9q1a6mtrWXdunWH3SvV1tbGK6+8wo9//OPD\n3nPVVVdx33338dprrw15vUM1ZsyYI06Ku91uSktLh72WoSotLT2i/46TtsHp9YP92+DIFowdzbqq\nqiouuugiPB4PABdeeGH4tUsvvbTfdbW1tXHuuecCcPnll/Pcc88Nef3R8nq9QPe5GMCRVzB6a3Xq\nVRin1w/2b4MjA8aOZt3RHkCXkZHR7/zxvqzu9Xrxer18+umn3HjjjXGtZagKCgoc9YXsy+n1g73b\n4MhDJDuadQsWLKCiooL29nb279/Ps88+e9T5s7OzycrKoqqqCoDVq1cPed1KJStHtmDsaNbNnj2b\nRYsWMWvWLE488URKSkoGbRE9/PDD4ZO8ZWVlQ163UslK4vFs6pKSEtN3wKl33nmHGTNmDHsth9q/\nfz+jR4/G7/dzzjnnsGrVKk4//fS41hSJRPjdqZFFRDYZY0oGm8+SFoyI/Cvwz4AB6oGrjDGOu0Fm\nyZIlvP3227S3t3PFFVc4IlyUSmQxB4yI5ALXAzONMQEReQr4NvCHWJc93B5//PF4l6BUUrHqJG8q\n4BGRVMAL7LZouUopB4s5YIwxLcCvgCZgD+AzxrwQ63KVUs4Xc8CIyFjgImAKMAHIEJHL+plviYjU\niEhNa2trrKtVSjmAFYdI5wO7jDGtxpggsA6Y13cmY8wqY0yJMaZk/PhBH6eilEoCVgRMEzBXRLzS\n3bW1FHjHguU6Tnl5OW+//Xb455/97Ge8+OKLA87/8ssvc8EFFwxHaUrFRcxXkYwxb4jIGmAz0AnU\nAqtiXe5gymtbWP58A7vbAkzI9rC0bBqLi3PtXu2AOjs7KS8v54ILLmDmzJkA3HHHHXGrR6lEYMlV\nJGPM7caY6caYfGPM5caYg1YsdyDltS3cuq6elrYABmhpC3DrunrKa1tiWm5jYyPTp0/niiuuoLCw\nkEsuuQS/388dd9zB7Nmzyc/PZ8mSJeH7ls477zx+8pOfcO6553LPPfewfv16li5dSlFRETt27ODK\nK69kzZo1AGzcuJF58+Yxa9Ys5syZE75JsdeBAwe4+uqrmT17NsXFxTz99NMAbNu2jTlz5lBUVERh\nYSHbt2+PaRuVGk6OvBdp+fMNBIKHj4USCHax/PmGmJfd0NDAkiVLqKurIzMzkwceeIDrrruOjRs3\nsnXrVgKBAM8880x4/t6hHG677TYWLVrE8uXL2bJlC1OnTg3P09HRwaWXXsq9997LW2+9xYsvvhi+\na7vXL37xCxYuXMjGjRvZsGEDS5cu5cCBA6xcuZIbbriBLVu2UFNTQ15eXszbqNRwcWTA7G4LRDU9\nGhMnTmT+/PkAXHbZZVRVVbFhwwbOPPNMCgoKeOmll9i2bVt4/v6GcuiroaGBE044gdmzZwOQmZlJ\naurhR6cvvPACd999N0VFRZx33nm0t7fT1NTEWWedxZ133sk999zDBx98cEQwKZXIHHmz44RsDy39\nhMmE7Ni/fH2HYBARfvCDH1BTU8PEiRNZtmzZYY8J6W8oh74iGdrBGMPatWuZNm3aYdNnzJjBmWee\nybPPPktZWRkPPvggCxcujGKLlIofR7ZglpZNw+N2HTbN43axtGzaAO+IXFNTU3iEuieeeCI8FGZO\nTg779+8Pn1Ppz5gxY444twIwffp0du/ezcaNG4HuQaI6OzsPm6esrIzf/va34fM7tbW1AOzcuZOT\nTjqJ66+/nkWLFlFXVxfzNio1XBwZMIuLc7nr4gJysz0IkJvt4a6LCyy5ijRjxgweeeQRCgsL+fTT\nT/n+97/PNddcQ0FBAYsXLw4f5vTn29/+NsuXL6e4uJgdO3aEp6elpfHkk0/yox/9iFmzZvHlL3/5\niIel/fSnPyUYDFJYWEh+fj4//elPAXjyySfJz8+nqKiId999l+9973sxb6NSw0WHazhEY2MjF1xw\nAVu3bo1rHdFKhN+dGlkiHa7BkS0YpZQzaMAcYvLkyY5rvSiVyDRglFK20YBRStlGA0YpZRsNGKWU\nbTRgDnHfffcxY8YMvvvd79qy/MbGRvLz821ZtlKJyJG3CgBQ9xRU3gG+ZsjKg9KfQeG3YlrkAw88\nwHPPPceUKVMsKlKpkc2ZLZi6p6DievB9CJjuvyuu754+RNdeey07d+5k0aJF/PrXv2bx4sUUFhYy\nd+7ccPf8ZcuW8atf/Sr8nvz8fBobG2lsbGTGjBlcc801nHbaaXzlK18hEOi+V2rTpk3MmjWLs846\ni/vvvz+mzVbKaZwZMJV3QLDPzY7BQPf0IVq5ciUTJkxgw4YNNDY2UlxcTF1dHXfeeWdE3fO3b9/O\nD3/4Q7Zt20Z2djZr164F4KqrruK+++4L39+k1EjizIDxNUc3PUpVVVVcfvnlACxcuJC9e/fi8/mO\n+p4pU6ZQVFQEwBlnnEFjYyM+n4+2tjbOPfdcgPAylRopnBkwWQMMujTQ9Cj1d3+WiJCamkooFApP\nO/SGxVGjRoX/7XK56OzsjGiYBqWSmTMDpvRn4O4z9ovb0z3dAueccw6rV68GugfmzsnJITMzk8mT\nJ7N582YANm/ezK5du466nOzsbLKysqiqqgIIL1OpkcKZV5F6rxZZfBWp17Jly7jqqqsoLCzE6/Xy\nyCOPAPDNb36TRx99lKKiImbPns2pp5466LIefvhhrr76arxeL2VlZZbUp5RT6HANSUB/d2q46XAN\nSqm404BRStkmoQImHodrTqe/M5XIEiZg0tPT2bt3r35homCMYe/evaSnp8e7FKX6lTBXkfLy8mhu\nbqa1tTXepThKenq6PoxNJSxLAkZEsoEHgXzAAFcbY6LqG+92u/UmQ6WSjFUtmHuBvxhjLhGRNMBr\n0XKVUg4Wc8CISCZwDnAlgDGmA+iIdblKKeez4iTvSUAr8LCI1IrIgyJyxPNURWSJiNSISI2eZ1Fq\nZLAiYFKB04HfG2OKgQPALX1nMsasMsaUGGNKxo8fb8FqlVKJzoqAaQaajTFv9Py8hu7AUUqNcDEH\njDHmI+BDEel98nwp8Hasy1VKOZ9VV5F+BKzuuYK0E7jKouUqpRzMkoAxxmwBBr2zUik1siTMrQJK\nqeSjAaOUso0GjFLKNhowSinbaMAopWyjAaOUso0GjFLKNhowSinbaMAopWyjAaOUso0GjFLKNhow\nSinbaMAopWyjAaOUso0GjFLKNhowSinbaMAopWyjAaOUso0GjFLKNhowSinbWPVUAdWP+vp6Kisr\n8fl8ZGVlUVpaSkFBQbzLithQ6i+vbWH58w3sbgswIdvD0rJpLC7Ojfj1oSqvbWHZ+m20BYIAjPW6\nuf3C05ia+qmjPwNw9n4kxphhX2lJSYmpqakZ9vUOp/r6eioqKggGg+FpbrebCy+80BE7x6H1BwIB\nQqEQbrebsrIyZs6c2e97/ly3h58/+zbtwa7Dpmd73Cwtmw5wxOvpbhcXzjqBqu17+cgX4PgsD9d9\n6WS+XngCf67bw+82vB+evuCUcVRt38seXwCXCF3GhP/uT2oKzE1rZiKf4PF4SElJcdRnAIm7H4nI\nJmPMoE8S0YCxyYoVK/D5fEdMz8rK4sYbb4xDRdHprX/r1q2sWbMm3uXErKioiMWLFwPO+Qwgcfej\nSANGD5Fs0t9OcbTpicbn87Gj8xieq92Fa/Q4MudeEn4tzZXColkTKJyYHZ5W92EbazY3x6PUowq1\n78dXtZqMjIzwNKd8BuD8/UgDxiZZWVkD/s/jBHvcufxfWwYHGt8i8/QLyTzjwsNe35btYdV1C4Hu\n8x/3r6sn84zieJR6VPvrXwQgc9rc8DSnfAbg/P3IsqtIIuISkVoRecaqZTpZaWkpbrf7sGlut5vS\n0tI4VRSdTZ157Ht/E3R14p0274jXW9oCTLnlWebf/RL/WbGNQJ/zLonC31CNK3M87+aczY7OYxz1\nGYDz9yMrWzA3AO8AmRYu07F6T8A58ex/eW0Lrf6u7i/n6HGkTZjW73yG7qBJVKGDfgKNtYwp/gYh\nSaW2ayJLL5zqiM+gl5P3I7AoYEQkD/gG8Avg36xYZjIoKChwzI7Qq7y2hVvX1RPqCNC+azOjZ5Uh\n4szuUoEdb/a0wBYAsD/kdtznAc7cj3pZteesAG4GQgPNICJLRKRGRGpaW1stWq2y2vLnGwgEuwjs\nqMF0duCdNj/eJR1VbraHsV53v68daKjGNfoYRuV2t8AmZHuGszSFBQEjIhcAHxtjNh1tPmPMKmNM\niTGmZPz48bGuVtlkd88hj7+hmpSMbEblzohzRUd34GAnn/mDR0wPdQRo37kJ77T5iKTgcbtYWtb/\noZ6yjxUtmPnAIhFpBP4ILBSRxyxYroqDCdkeQsF2Ajs34j11HpLiGpb1jvW68bqj3x17e+72FW6B\nnTqP3GwPd11cYEmPYRWdmAPGGHOrMSbPGDMZ+DbwkjHmspgrU3GxtGwaoaYtmODBYT08+swfZGzG\nKMuW53/vVVK82WRPKbDsdgQVPWeevVO2WVycy9QDW0nxZJI+MT/q97tEhrTeFPni8CxWoWA7gR0b\n8Z56Fu1d3eeVVHxYGjDGmJeNMRdYuUw1vNrb26mtqiQnf8GQDo8yPUO7MBky3Ze9rdC+azMm2B5u\ngVkVXCp62oJRh/nrX//Kvn37GD19aIdH/Z1wHW4HGqq7W2CTui/t6tWj+NGAUYe57b/+m5T00XQe\nd1q8SxkS0xkk8P6beE+ZG26B6dWj+NGAUWH/++Yutr72Ep6T5yIuZ96mFmisxXQEwodHY71uPcEb\nRxowKuz23z+JOXig33uPnMLfUE3KqAzSTywE4PYLndkSSxYaMCrsw80bkDQPnsmJd1d0JExXkMD2\n1/GcMhdxufG4U7T1EmcaMAqAYDBI+/bX8J58JpLaf9f7RNf+QR2hQ1pgd11cGOeKlAaMAuCVV16h\nM7CPMTMWxLuUIfM3VIdbYO4UtPWSADRgFABr164lIyODzJPPiHcpQ2JCXfi3v47n5DlIahrL/7Eo\n3iUpNGAU0NXVxbp16yiat5CDDP3waFFKFVVp17Nz1D9RlXY9i1KqLKzy6Nqb6gkFPiej5+qRtl4S\ngzOvRSpLVVVV8fHHH+PJHvo5i0UpVdztfhCvdACQJ59wr/sBlvEoFV1zKU3ZwgT5hN0mh192fov1\nIWsPxfwN1Yh7FOlTTh/y7QrKehowijVr1pCWnk4or2jITdqbU58Kh0svETiG/XzP9SK93/k8+YS7\n3Q9CEMtCpvvw6DU8J80mxZ0+4GNM1PDTQ6QRLhQK8fiT/4t7UjEpaYN3qU9N6b91MEE+GfA9fRsU\nXung5tSnoqrzaA62vEPoQFu4c522YBKHBswI9/rrr/Np69/xRDg0Q1fI4HEffhOkx+3iI8mJar0T\nZG9U8x+Nv6EaSU3DM7X7MT3agkkcGjAj3Jo1a8CVimfqnIjmN8BdFxeQm+1BIDyYU8vpNxMwaRGv\nd7cZN7SC+9ZjQvgbqkmfcnq4BZarNzcmDD0HM4IZY1i7di2eKaeTMsob0XsuclUx9+kb+D/Tysfp\n4/lw5lL4cDsTNy8nnQ6MOfKQ6Mj1QmXImsvIHbsb6Nr/afjwSIfGTCwaMCNYTU0NTU1NTPqHb0U0\n/6KUKu5KfRAvHSBwPK2M3XQrgpAmnRDhqQ8RKE3Zwu0x1N7rQEM1uFLxnjyH3GyPjl6XYPQQaQRb\ns2YNqamp/OKGKyOav78rRaOkqztcomTFORhjDP6GajyTi8kYnanhkoA0YEYoYwxr1qzh/PPP57Jz\nTyPbM3gHu6NdKYpWGxmDzzSIjo+20/V5K95p8wkEu3RozASkATNCbdmyhZ07dzJ1Tinz735pwNH5\nD7XbRHel6GjGsp/3R13Gf6Y+NORl+BtehRQXnpPPBHRozESkATNCrV27Fklxsf6zCRE//vWXnd/C\nH8WVoqMRgVQJ8T3Xi0MKmd7Do/QTZ+HyjAF0aMxEpAEzAhljePixJxg1MZ8Ub1bE71sfWsBPuq5h\nNzmEjNBlYu/QJgLfdb0U9fuCH++is22PXj1KcBowI9C2bdvY/cHOIT33KBQydPdjM1j1HADXwE8c\nHtCBhmqQFLynzMUlog9WS1AaMCNMeW0LX7txOSB4Tz0rqvf23tCYK5+QIuCysEd+NHdehw+PJuUz\nOusYfv2tWRouCUoDZgQpr23h1nX1/P2tlxk18TRcGWOjen9/l6mtIEJU9yZ17W2i89NmJhR9SVsu\nCU472o0gy59v4POPGgl+0sTY8/8l6vdbeZn6yGVH3i8mpfFNRITqlbdw/PHH21aTil3MLRgRmSgi\nG0TkHRHZJiI3WFGYsl5LW6D70i4M6fAoNMDuYsW9hdHcm/TRWy+zYMECDRcHsOIQqRP4sTFmBjAX\n+KGIzLRgucpiLhEONFQzKncGqWMi79PSe+4lVfo/GRvr6AgdJpVfdkZ2u0JwbzPBTz7gkksuiW2l\naljEHDDGmD3GmM09/94HvAPoQXECav90N8GPd+I9NbrnHtl17qXXPpMe8eBT/ve6W2AXX3yxbfUo\n61h6kldEJgPFwBtWLldZw930JkDUl6ftPPcCMFYORDyvv6GaMZNmkpeXZ2NFyiqWBYyIjAbWAjca\nYz7v5/UlIlIjIjWtra1WrVZFwf3hm6RPOJXUrGOjep+Vtwj0v/zIzr8EP9tDx993kDo1uvNHKn4s\nCRgRcdMdLquNMev6m8cYs8oYU2KMKRk/frwVq1URKq9t4Yx/f4L3t73F2NPOJtvjRoj83EmktwgY\nE/0JX79Ji/j8i/+9agCmzF4Y3UpU3MR8mVpEBPgf4B1jzH/FXpKyUrjvy8ZKAGTKXHyBIOHOuBFY\nH1oAQbjX/cBRQ0kkuoAxBm4J/nPk518aXiXt+FP42lkFka9ExZUVLZj5wOXAQhHZ0vPn6xYsV1lg\n+fMNBIJd+N97FfexJ+Eee8KQO/hH8r7+Amig0DnAqIjDpdP3MR173sM7bT7P1u2J6D0q/mJuwRhj\nqoh4LDM13Ha3Bejc9wkHW94h++zLh7ycm1OfYoAHCgxqoFZPahT3IPVePfJOm8dn/sGHllCJQW8V\nSHITsj3433sNiP7q0WHLseFK0iiCEQ10Bd1Xj9zHTsE9doLldSj7aMAkuaVl07q/nDmTcI/r/9Ju\n73OEjtZAsetK0sHOwVsxvS2w3oCMNJRU/GnAJLm2va0c/HDbUVsv3zlzIo13f4PfXDrwSP9WDjbV\nKwQEgl2DztfbAsuYNh93irBs0WmW1qHsozc7Jrn/WPEQYI4aMKvfaGL1601MyPbgdafgDx7Zqoj0\nSlKkjIHHus6PaN6uHa/jzpnE5Kmn6sDeDqMtmCRWXtvC3q1/I/WYPNw5Jw44nzHdV4ha2gL9hksv\nKx9YfxAXt3dePeh8xt/GvsY6bv3BlVTfslDDxWE0YJLYnevepL2pHu+0eYhFz2v+jNExLyNohJuD\nkQ0Xsb/hVYwxenOjQ2nAJLEdNS+DCZERw9WjvpYFv0eHie3Ieh8ZUXWuG3PcJPLz82Nap4oPDZgk\nVV7bQqChmtTs43Efe5Jly10fWsBNwSU0h3KGPA5MNpHd3Njl99HeVIecNNeyFpgaXhowSequ8hoC\nH2zBO22+5V/O9aEFLOi4jxuCPxjSlaVIb270b38DTCjq4SVU4tCrSElqZ80rEOqK+svpdgnBrsia\nJr1Xlm5OfYpc+WTA2wQOnR7VzY0N1aRmHUfacVMjml8lHm3BJKH/KK/nwHvVuMaMJ+2EUyN+X262\nh0tnTwx3vIvEoa2ZvudmOkwqj3adT3Oo+zlKzaGciG9u7GrfT/sHb4VbYOW1LRHXpBKHtmCSzH+U\n1/PoK+8S2LWZMcXfiPjwaEVPJ7tb19XTNYSTK4e2ZibIXnabcfyy81usDy3g9qiXBoH334BQZ7j/\nzvLnG/QStQNpwCSZ1a83EdjxJnR14p0Web+VG5/cQopAKIYBvNeHFrC+w5q+Mv6Gw1tg+txpZ9JD\npCRj6H7qoWv0MYzKje5RqrGEi5VCB/0Edm0+rP+OPnfamTRgkkyoI0D7zk14T52HiDM/3i9aYF/0\n39HnTjuTM/dANSDTVIvp7IhpaIZ48ze82tMCmw7AZXMn6fkXh9JzMEnmxM/r2ePNZlSecx5NlSKQ\nme7GFwgyxtVJ084aRhd+mbyxGXpzo8NpwCQRv9/PpqpKsmeeh6S44l1OxEIG2gJBcrM9zEvZTn1n\nB+t/fRPnnXdevEtTMdKASSJ3/P5x2gN+Mk92Zs/XlrYAK9evJnPsOM4+++x4l6MsoOdgksj/PPZH\nUjyZpE+0/sbA3GwPuTZfyQkFD3Lg/Y2kTJmDy+WcFpgamAZMkjh48CB7334V7ylzEZf1DdMvTR/P\n0rJpeNz2ffHbd23GBNtJO3me9txNEhowSeKOlX/EdAT6vfcoI81FbrYHobslctncSVEvf8O7rSwu\nzuWbZ9h3wtXfUE1K+hjSJxWw/PkG29ajho+eg0kC5bUt/O6hx5BRGaRPnnXE625XCtW3HP40xA3v\nttISRe/Y3W0BymtbWLvJnpaF6Qzif//N7nuPXKnaczdJaAsmCdzz7Fb2N7yG95QzEdeRI+77Akc+\nR2hp2TTcrshvaszyuPnXJ7dENEj3UAQ+2ILp8IcHx9Keu8lBAyYJ7Kp7g9DBAwN2ruvvy7q4OJeM\ntMgasO4U4fP24JCfCBkJ/7vV4RaYx+3SnrtJQg+RkoBpfANJ8+CZXNzv6y1tAabe+me6jCE328OX\npo9nw7uttPXTsulLBII236RkuoIEtr+G9+Q5pKS6ueviAu1clyQsCRgR+SpwL+ACHjTG3G3Fcu1U\nX19PZWUlPp+PrKwsSktLKShw1kPV6+vreeGFF/BtfQXP1NlI6sCjy/UOwdDSFuCx15siXsdQh8WM\nRvsHdeEW2L+ddYyjwiUZ9iM7xXyIJCIu4H7ga8BM4DsiktD91Ovr66moqMDn8wHg8/moqKigvr4+\nzpVFrncb6urqOBg4YOnA3sPN/96rSJqHwql57Nu6wTGfQzLsR3azogUzB3jfGLMTQET+CFwEvG3B\nsm1RWVlJMHj44UEwGKSystIx//v0bsOjjz4KwP6acvZt+Uucqxqag7vf4cRTTuPsjI8IBnHM55AM\n+5HdrAiYXODDQ35uBs7sO5OILAGWAEyaFH0/DCv1/o8T6fRE1LfWDBNgX0fI1hOxVkvBMFo6OPb4\nYymd+8Vja53yOSTDfmQ3KwKmv2udR+znxphVwCqAkpKSuH4PsrKy+t0JsrKy4lDN0PRuw7Jly8LT\ndnQewxvBSRy06dy9x+3irou7/2de/nxDVP1oegnwm0uLWFycy4oVKxz9OSTDfmQ3Ky5TNwMTD/k5\nD9htwXJtU1paitt9eH8Rt9tNaWlpnCqKXn/bMN2zjzXfPYkVlxaF7xvqHcC79+/enrzZniP7y7hd\nwmVzJ4V7/Y71usn2uMM9gHuv7iwuzqX6loWsuLQId0rkfWkE+O4hY7s4/XNwev3DQUyMlwlEJBV4\nDygFWoCNwD8ZY7YN9J6SkhJTU1MT03pjlQxn/2PdhvLaFpY/38DutgATsj1DGnul7zJ6L4HvbguQ\n7XVjTHdHv4GW7/TPwen1D5WIbDLGlAw6X6wB07OyrwMr6L5M/ZAx5hdHmz8RAkYpNXSRBowlB+vG\nmD8Df7ZiWUqp5KG3CiilbKMBo5SyjQaMUso2GjBKKdtowCilbKMBo5SyjQaMUso2GjBKKdtowCil\nbKMBo5SyjQaMUso2GjBKKdtowCilbKMBo5SyjQaMUso2GjBKKdtowCilbKMBo5SyjQaMUso2GjBK\nKdtowCilbKMBo5SyjQaMUso2GjBKKdtowCilbBNTwIjIchF5V0TqRORPIpJtVWFKKeeLtQXzVyDf\nGFMIvAfcGntJSqlkEVPAGGNeMMZ09vz4OpAXe0lKqWRh5TmYq4HnBnpRRJaISI2I1LS2tlq4WqVU\nokodbAYReRE4vp+XbjPGPN0zz21AJ7B6oOUYY1YBqwBKSkrMkKpVSjnKoAFjjDn/aK+LyBXABUCp\nMUaDQykVNmjAHI2IfBX4d+BcY4zfmpKUUski1nMwvwPGAH8VkS0istKCmpRSSSKmFowx5mSrClFK\nJR/tyaue1HDiAAAEwUlEQVSUso0GjFLKNhowSinbaMAopWyjAaOUso0GjFLKNhowSinbaMAopWyj\nAaOUso0GjFLKNhowSinbaMAopWyjAaOUso0GjFLKNhowSinbaMAopWyjAaOUso0GjFLKNhowSinb\naMAopWyjAaOUso0GjFLKNhowSinbaMAopWyjAaOUso0GjFLKNpYEjIjcJCJGRHKsWJ5SKjnEHDAi\nMhH4MtAUezlKqWRiRQvmN8DNgLFgWUqpJBJTwIjIIqDFGPNWBPMuEZEaEalpbW2NZbVKKYdIHWwG\nEXkROL6fl24DfgJ8JZIVGWNWAasASkpKtLWj1AgwaMAYY87vb7qIFABTgLdEBCAP2Cwic4wxH8Va\nWH19PZWVlfh8PrKysigtLaWgoCDWxQ4rp2+D0+uH5NgGJxs0YAZijKkHju39WUQagRJjzCexFlVf\nX09FRQXBYBAAn89HRUUFgGN2Dqdvg9Prh+TYBqdLyH4wlZWV4Z2iVzAYpLKyMk4VRc/p2+D0+iE5\ntsHphtyC6csYM9mqZfl8vqimJyKnb4PT64fk2AanS8gWTFZWVlTTE5HTt8Hp9UNybIPTJWTAlJaW\n4na7D5vmdrspLS2NU0XRc/o2OL1+SI5tcDrLDpGs1HsCzsln/52+DU6vH5JjG5xOjBn+LiklJSWm\npqZm2NerlLKGiGwyxpQMNl9CHiIppZKDBoxSyjYaMEop22jAKKVsowGjlLKNBoxSyjYaMEop28Sl\nH4yItAIfDPuK+5cDxHwHeJzpNsSf0+uH6LbhRGPM+MFmikvAJBIRqYmkw1Ai022IP6fXD/Zsgx4i\nKaVsowGjlLKNBkzPOMEOp9sQf06vH2zYhhF/DkYpZR9twSilbKMBo5SyjQYMICLLRKRFRLb0/Pl6\nvGuKlIh8VUQaROR9Ebkl3vVES0QaRaS+5/fuiEGCROQhEflYRLYeMu0YEfmriGzv+XtsPGsczADb\nYPn3QAPmC78xxhT1/PlzvIuJhIi4gPuBrwEzge+IyMz4VjUkX+r5vTulH8kfgK/2mXYLUGmMOQWo\n7Pk5kf2BI7cBLP4eaMA42xzgfWPMTmNMB/BH4KI415T0jDF/Az7tM/ki4JGefz8CLB7WoqI0wDZY\nTgPmC9eJSF1P0zGhm7eHyAU+POTn5p5pTmKAF0Rkk4gsiXcxMTjOGLMHoOfvYweZP1FZ+j0YMQEj\nIi+KyNZ+/lwE/B6YChQBe4Bfx7XYyEk/05zW72C+MeZ0ug/zfigi58S7oBHM8u9BQj5VwA4DPWO7\nLxH5b+AZm8uxSjMw8ZCf84DdcaplSIwxu3v+/lhE/kT3Yd/f4lvVkPxdRE4wxuwRkROAj+NdULSM\nMX/v/bdV34MR04I5mp4dotc/AFsHmjfBbAROEZEpIpIGfBtYH+eaIiYiGSIypvffwFdwzu++r/XA\nFT3/vgJ4Oo61DIkd34MR04IZxC9FpIjuw4tG4F/iW05kjDGdInId8DzgAh4yxmyLc1nROA74k4hA\n9774uDHmL/EtaXAi8gRwHpAjIs3A7cDdwFMi8v+AJuAf41fh4AbYhvOs/h7orQJKKdvoIZJSyjYa\nMEop22jAKKVsowGjlLKNBoxSyjYaMEop22jAKKVs8/8BBx7np9jFbjEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6c9c0e77b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=plt.figaspect(1))\n",
    "ax.scatter(grid[:, 0], grid[:, 1], color='#808080', label='grid')\n",
    "ax.scatter(pos[:, 0], pos[:, 1], label='particles')\n",
    "ax.scatter(pos[neighbors, 0], pos[neighbors, 1], label='found')\n",
    "\n",
    "ax.plot([0, a[0], a[0] + b[0], b[0], 0], [0, 0, b[1], b[1], 0], color='black', zorder=4)\n",
    "ax.legend()\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:mda]",
   "language": "python",
   "name": "conda-env-mda-py"
  },
  "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.2"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}