{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyzing GROMACS data with freud and MDTraj: Computing an RDF for Water\n",
    "\n",
    "In this notebook, we demonstrate how `freud` could be used to compute the RDF of the output of an atomistic simulation, namely the simulation of TIP4P water.\n",
    "In the process, we show how the subsetting functionality of such tools can be leveraged to feed data into `freud`.\n",
    "We use this example to also demonstrate how this functionality can be replicated with pure NumPy and explain why this usage pattern is sufficient for common use-cases of `freud`.\n",
    "The simulation data is read with [MDTraj](http://mdtraj.org/) and the results are compared for the same RDF calculation with `freud` and MDTraj."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating water\n",
    "\n",
    "To run this notebook, we have generated data of a simulation of TIP4P using [GROMACS](http://www.gromacs.org/).\n",
    "All of the scripts used to generate this data are provided in this repository, and for convenience the final output files are also saved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import freud\n",
    "import mdtraj\n",
    "import numpy as np\n",
    "\n",
    "traj = mdtraj.load_xtc(\"output/prd.xtc\", top=\"output/prd.gro\")\n",
    "bins = 300\n",
    "r_max = 1\n",
    "r_min = 0.01\n",
    "\n",
    "# Expression selection, a common feature of analysis tools for\n",
    "# atomistic systems, can be used to identify all oxygen atoms\n",
    "oxygen_pairs = traj.top.select_pairs(\"name O\", \"name O\")\n",
    "\n",
    "mdtraj_rdf = mdtraj.compute_rdf(traj, oxygen_pairs, (r_min, r_max), n_bins=bins)\n",
    "\n",
    "# We can directly use the above selection in freud.\n",
    "oxygen_indices = traj.top.select(\"name O\")\n",
    "\n",
    "# Alternatively, we can subset directly using Python logic. Such\n",
    "# selectors require the user to define the nature of the selection,\n",
    "# but can be more precisely tailored to a specific system.\n",
    "oxygen_indices = [atom.index for atom in traj.top.atoms if atom.name == \"O\"]\n",
    "\n",
    "freud_rdf = freud.density.RDF(bins=bins, r_min=r_min, r_max=r_max)\n",
    "for system in zip(np.asarray(traj.unitcell_vectors), traj.xyz[:, oxygen_indices, :]):\n",
    "    freud_rdf.compute(system, reset=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fb10d1898e0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEXCAYAAACpuuMDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3RU55nv+e9TqpJUSEISkkCypEKYiy3uF4FlcDA2lzjG2E7H03EnISdZx+24u5N1TibxdLpndTLT56yenDl9OieOO3GYbo/bnrSdOL7EODYBO8Y2xrIBg7kJgwAhJCQQEqVr6VJVz/xRhawSwpRQqUqX57NWLVS7dpWejfH+1bvfd7+vqCrGGGPMZY5EF2CMMWZ0sWAwxhgTwYLBGGNMBAsGY4wxESwYjDHGRLBgMMYYE8GCwRhjTAQLBmOGSESqRcQnIu0i0iAiT4lIevi1p0SkR0Tawo/DIvJ/iUhmv/d/Q0QC4fdffjyeuCMyJpIFgzHXZ5OqpgOLgSXA3/R77f9W1QwgD/gmUA68JyJp/fZ5X1XT+z2+HbfKjbkGCwZjhkFVG4A/EAqIga91qeoe4F4gh1BIGDPqWTAYMwwiUgR8Aai62j6q2gbsAD4Xr7qMGQ4LBmOuz8si0gacBS4AP7rG/ueAKf2el4uIt9+jfKQKNWaoLBiMuT73h/sR1gA3A7nX2L8QaO73vEJVs/o9KkaoTmOGzILBmGFQ1beBp4B/vNo+4RFL64B341SWMcPiTHQBxowD/xOoFpGIDmgRSQHmA/8NuAT8vwmozZghsxaDMcOkqo3A08DfhTf9b+H+h+bw9n3ASlXtSFCJxgyJ2EI9xhhj+rMWgzHGmAgWDMYYYyJYMBhjjIlgwWCMMSbCmB+umpubqyUlJYkuwxhjxpR9+/ZdVNW8wV4b88FQUlLC3r17E12GMcaMKSJy5mqv2aUkY4wxESwYjDHGRLBgMMYYE2HM9zEYY0x/vb291NbW0tXVlehSRoXU1FSKiopwuVxRv8eCwRgzrtTW1pKRkUFJSQkikuhyEkpVaWpqora2lhkzZkT9PruUZK5PTQVcqo7cdqk6tN2YBOrq6iInJ2fChwKAiJCTkzPk1pO1GMz1ycjnfMVzvOq/hcquKZSmNnOP8wOmlT+Y6MqMsVDo53r+LqzFYK5LZVc2Wy4uwlO/nRUcxlO/nS0XF1HZlZ3o0owxw2TBYK7LtsPnCWR68GYvoLh1P97sBQQyPWw7fD7RpRmTcI899hilpaV89atfjenn7ty5k3vuuSemnzkYu5Rkrkud18dNKU3ktx+ldvJS8tuP4k3J5xNvTqJLM2ZIKutb2Hb4PHVeH4VZbu6aP43SgsxhfebPf/5zXn/99YgOX7/fj9M5Nk651mIw16U0tZlFZ57mXMZ8arPKOJ67lpKGHSyXSuuANmNGZX0LW945TYuvl4LMVFp8vWx55zSV9S3X/ZmPPPIIp06d4t577yUzM5OHH36YDRs28PWvf51AIMCjjz7K8uXLWbhwIb/85S+BK1sC3/72t3nqqacA2LZtGzfffDO33XYbL7744rCON1pjI77MqLO2KMDWxjvwNOyjqiGZT7pzWBQs4jspb0HG3yS6PGOisu3weTLdLjLdoTH+l//cdvj8dbcannjiCbZt28Zbb73F448/ztatW9m1axdut5stW7aQmZnJnj176O7uZtWqVWzYsOGqn9XV1cWf//mf88c//pFZs2bx5S9/+bpqGiprMZjrUrJkHTeWrePF7uUsat3J55xHWe2u5hnH/dYBbcaMOq+PjNTI78cZqU7qvL6Y/Y57770Xt9sNwPbt23n66adZvHgxt9xyC01NTZw4ceKq7z127BgzZsxg9uzZiAhf+9rXYlbXZ7EWgxm6mgrIyOf4+R5mzJ5HdrefsuadtKQW93VAD/carTHxUJjlpsXX29dSAGjr8lOY5Y7Z70hLS+v7WVX52c9+xuc///mIfXbt2kUwGOx73v++g0QMvbUWgxm68D0MBw99zOkTh3FVv4W/u5O0ngsUcj6m37aMGUl3zZ9Gi6+XFl8vQdW+n++aP21Eft/nP/95fvGLX9Db2wvA8ePH6ejoYPr06Rw9epTu7m5aWlp48803Abj55ps5ffo0J0+eBODZZ58dkboGshaDGbLKrmx+e3ERG3t+S7q20y7p/HPvvSyYOplFDTtwFFz9mqkxo0lpQSYPr54RMSrpy8uLRqzF+9BDD1FdXc3SpUtRVfLy8nj55ZcpLi7mT//0T1m4cCGzZ89myZIlQGieoy1btrBx40Zyc3O57bbbOHz48IjU1p+o6oj/kpFUVlamtlBPfP1kx3FafL0sa/g1/oun2Jd6KwddCxERyjLb2DzfRcmSdYku00xQlZWVlJaWJrqMUWWwvxMR2aeqZYPtH7dLSSKSKiIfisjHInJERP7PQfYREXlMRKpE5KCILI1XfSZ6dV4fhZwnyxWkqWgdN+tJpvTU0xMI8sC6lRYKxoxx8byU1A3cqartIuICdonI66raf9D7F4DZ4cctwC/Cf5pRpDS1GU/9Do4XbKAttRDJns5X67dTU7DBOp2NGQfi1mLQkPbwU1f4MfA61n3A0+F9K4AsESmIV40mOmuLAlS4V3NWpxJU5axOpcK9mrVFgUSXZoyJgbiOShKRJBE5AFwAdqjqBwN2KQTO9nteG95mRpGSJet4YN1KMt0u6lu6yHS77BKSMeNIXEclqWoAWCwiWcBLIjJfVft3sQ82YPeK3nEReRh4GMDj8YxIreazlRZk2mUjY8aphNzHoKpeYCdw14CXaoHifs+LgHODvH+LqpapalleXt6I1WmMMRNRPEcl5YVbCoiIG1gHHBuw2yvA18Ojk8qBFlWtj1eN5jrYSm7GxERJSQkXL17E6/Xy85///Lo+4+6778br9Q67lnheSioA/k1EkggF0m9U9VUReQRAVZ8AXgPuBqqATuCbcazPXA9byc2MZeHpXcgu+XTbpWpoawBPeUJKuhwMf/mXf3nFa4FAgKSkpKu+97XXXotJDfEclXRQVZeo6kJVna+qfx/e/kQ4FC6PXPorVZ2pqgtU1e5cG+VsJTczpmXkw5GXP231XqoOPc/IH9bHVldXc/PNN/PQQw8xf/58vvrVr/LGG2+watUqZs+ezYcffkhTUxMbNmxgyZIlfOtb3+LyzcY/+MEPOHnyJIsXL+bRRx9l586d3HHHHXzlK19hwYIFANx///0sW7aMefPmsWXLlr7fe7nVMVw2JYYZlr6V3FIXUNz6EbXZSwmk2ER6ZozILoF594fCoHAp1H0Uet6/BXGdqqqqeP7559myZQvLly/n3//939m1axevvPIK//AP/4DH4+G2227jhz/8Ib///e/7TvA//vGPOXz4MAcOHABCazV8+OGHHD58uG/hnyeffJIpU6bg8/lYvnw5X/rSl8jJid0iWRYMZlhsJTcz5mWXhEKh+j0oWRWTUACYMWNG3zf8efPmsXbtWkSEBQsWUF1dTXV1dd/COxs3biQ7++qt7BUrVkSsBvfYY4/x0ksvAXD27FlOnDhhwWBGj/53QbemFtKaWkBJ/XabSM+MHZeqQy2FklWhP7M8MQmHlJSUvp8dDkffc4fD0bfMZ7RTavefunvnzp288cYbvP/++0yaNIk1a9ZETNMdCzbtthkWuwvajGmX+xTm3Q8zVn96WWngSLsRsHr1an71q18B8Prrr3Pp0iUAMjIyaGtru+r7WlpayM7OZtKkSRw7doyKitiPALRgMMNid0GbMa2tIbJP4XKfQ1vDiP/qH/3oR7zzzjssXbqU7du3992sm5OTw6pVq5g/fz6PPvroFe+766678Pv9LFy4kL/7u7+jvDxy9FQsFvaxabeNMePKRJ12OxAIMHXqVBoaGnC5XBGvjdppt40xxoycefPm8dBDD10RCtfDOp+NMWYcOHZs4EQS189aDMaYcWesXyKPpev5u7BgMMaMK6mpqTQ1NVk4EAqFpqYmUlNTh/Q+u5RkhmYUzi1jTH9FRUXU1tbS2NiY6FJGhdTUVIqKiob0HgsGMzQ2aZ4Z5VwuV8Rdwmbo7FKSGRKbNM+Y8c+CwQxJ36R52Qsobt2PN3sBgczQpHnGmPHBgsEMSZ3XRyHnIybNK+Q8dV5fokszxsSI9TGYIbFJ84wZ/6zFYIbEJs0zZvyzYDBDMtikeZvnuygpmRW5o637bMyYZZeSzJCVFmRGrs52KdmGsBozjliLwQybDWE1ZnyxYDDDZkNYjRlfLBjMsNkQVmPGl7gFg4gUi8hbIlIpIkdE5D8Nss8aEWkRkQPhxw/jVZ+5fqWpzZQ07OB47lpqs8o4nruWkoYdlKY2J7o0Y8x1iGfnsx/4nqp+JCIZwD4R2aGqRwfs966q3hPHuswwrS0K8Ix3NQGdSkZ4COs592o22xBWY8akuAWDqtYD9eGf20SkEigEBgaDGWNKlqzjgfwWth0OXT4qzHJz1/KVlPQfuWSMGTMSMlxVREqAJcAHg7x8q4h8DJwDvq+qR+JYmrlOVwxhNcaMWXEPBhFJB14A/rOqtg54+SNguqq2i8jdwMvA7EE+42HgYQCPxzPCFRtjzMQS11FJIuIiFAq/UtUXB76uqq2q2h7++TXAJSK5g+y3RVXLVLUsLy9vxOs2xpiJJJ6jkgT4V6BSVf/pKvvkh/dDRFaE62uKV43GGGPieylpFbAZOCQiB8Lb/hbwAKjqE8ADwF+IiB/wAQ+qLdxqjDFxFc9RSbsAucY+jwOPx6ciY4wxg7E7n40xxkSwYDDGGBPBgsEYY0wECwZjjDERLBiMMcZEsGAwsVFTEVrOsz9b3tOYMcmW9jSxkZFvy3saM05Yi8HEhC3vacz4YcFgYsKW9zRm/LBgMDFhy3saM35YH4OJidLUZjz1OzhesIHW1EJaUwsoqd+Oo2BDokszxgyRtRhMTKwtClDhXs1ZnUowvLxnhXs1a215T2PGHGsxmJiw5T2NGT8sGEzM2PKexowPdinJGGNMBGsxmOjVVFDlS2NrTUrf5aJNnm5muTvAU57o6owxMWItBhO1Kl8alW8+g3jPUJCZinjPUPnmM1T50hJdmjEmhiwYTNS21qRwMm89y9rfxtOyj2Xtb3Mybz1ba1ISXZoxJoYsGEzU6rw+ApkeGtLnUtT6EQ3pcwlkeuwmNmPGGQsGE7XCLDdJLTURdzcntdRQmOVOdGnGmBiyzmcTtU2ebipP7GBf3vpQS4E8ZjbuoHTh5kSXZoyJIQsGE7VZ7g5Yu5mqmhTqvT4Ks6ZTunBzaLsxZtyIWzCISDHwNJAPBIEtqvrTAfsI8FPgbqAT+IaqfhSvGs01eMqZBXz3pkQXYowZSfFsMfiB76nqRyKSAewTkR2qerTfPl8AZocftwC/CP9pjDEmTuLW+ayq9Ze//atqG1AJFA7Y7T7gaQ2pALJEpCBeNRpjjEnQqCQRKQGWAB8MeKkQONvveS1Xhgci8rCI7BWRvY2NjSNVphkqW/fZmHEh7p3PIpIOvAD8Z1VtHfjyIG/RKzaobgG2AJSVlV3xukkQW/fZmHEhri0GEXERCoVfqeqLg+xSCxT3e14EnItHbWb4bN1nY8aHuAVDeMTRvwKVqvpPV9ntFeDrElIOtKhqfbxqNMNj6z4bMz7E81LSKmAzcEhEDoS3/S3gAVDVJ4DXCA1VrSI0XPWbcazPDFOd18dNKU0Rd0Z7U/L5xJuT6NKMMUMQt2BQ1V0M3ofQfx8F/io+FZlYs3WfjRkfbK4kEzO27rMx44NNiWFixtZ9NmZ8sGAwMWXrPhsz9tmlJGOMMREsGIwxxkSwYDDGGBPBgsEYY0wECwZjjDERLBiMMcZEsGAwxhgTwYLBGGNMBAsGY4wxEYZ857OIpAFdqmoT4Jgr1VRQ5Utja01K37QYmzzdzHJ3gKc80dUZY6JwzRaDiDhE5Csi8nsRuQAcA+pF5IiI/HcRmT3yZZqxosqXRuWbzyDeMxRkpiLeM1S++QxVvrREl2aMiVI0l5LeAmYCfwPkq2qxqk4FPgdUAD8Wka+NYI1mDNlak8LJvPUsa38bT8s+lrW/zcm89WytSUl0acaYKEVzKWkdEAB+oKoHL29U1WZCy3S+EF6y0xjqvD4KMj00MJei1o+onbyUQKaHeq8v0aUZY6J0zRaDqvaqapBQQFx1n5hWZcaswiw3SS01Eau4JbXUUJjlTnRpxpgoDWVU0n4R+ZGI2Egmc1WbPN3MbNzBvvTbqclcxr7025nZuINNnu5El2aMidJQRiUVAwuAvxCRD4CDwEFVfX5EKjNj0ix3B6zdTFVNCvVeH4VZ0ylduDm03RgzJkQdDKr6pwAikgLMIxQStwAWDOZTnnJmAd+9KdGFGGOu1zWDQUREVfXyc1XtBj4KPwbdxxhjzNgV1XBVEfmOiHj6bxSRZBG5U0T+DfgPI1OeMcaYeIsmGO4iNFz1WRGpF5GjInIaOAH8GfATVX3qWh8iIk+KyAUROXyV19eISIuIHAg/fjiE4zDGGBMj17yUpKpdwM+Bn4fvV8gFfKrqHeLvegp4HHj6M/Z5V1XvGeLnGmOMiaGoO59F5ARwCPgYOCAiB1T1TLTvV9V3RKRkyBUaY4yJq6EMV/0lcCPQBHwB+FX4ktJLwH+J0U1ut4rIx8A54PuqemSwnUTkYeBhAI/HM9guJtZscjxjJoyhBMPXVHXx5Sci8gTwTaAV+CfgO8Os5SNguqq2i8jdwMvAoBP0qeoWYAtAWVmZjYaKg77J8fLWU5DpCU2Od2IHrN3MrIE7W4gYM6YN5S7mFhFZePmJqh4AylX1H4FVwy1EVVtVtT3882uAS0Ryh/u5JjaGMjmezbBqzNg2lBbDtwhdPjoAHABuAoLh15KHW4iI5APnVVVFZAWh0Goa7uea2BjK5Hhba1KQcIg0MJf89qPsy1tPVU2K3fhmzBgwlDufj4VP2H8CLASqgB+FF+557lrvF5FngTVArojUAj8CXOHPfgJ4gNB0G37ABzxoN82NHoVZbsR7JmJyvDryKMyafsW+NsOqMWPbkFZwC6/a9jxXToPxX6N4759d4/XHCQ1nNaPQJk83lSd2sC9vPYFMD3XkMbNxB6ULN1+x71BCxBgz+gx5aU8zMQ1lcryhhIgxZvSxYDDRGcLkeDbDqjFjmwWDiT2bYdWYMc0W3THGGBPBgsEYY0wECwZjjDERLBjMyKipgEvVkdsuVYe2G2NGNet8NiMjI5/zFc/xqv8WKrumUJrazD3OD5hW/mCiKzPGXIO1GMyIqOzKZsvFRXjqt7OCw3jqt7Pl4iIqu7ITXZox5hosGMyI2Hb4PIFMD97sBRS37sebvYBApodth88nujRjzDVYMJgRUef1Ucj5iGkxCjlPnc2XZMyoZ30MZkSUpjbjqd/B8YINtKYW0ppaQEn9dhwFGxJdmjHmGqzFYEbE2qIAFe7VnNWpBFU5q1OpcK9mbVEg0aUZY67BWgxmRJQsWccD+S1sO3y+bxW3u5avpKQgM9GlRaqpgLZ6AKqcs9hak4LvwkkWOk6zuDiTIs9MW3XOTDgWDGbElBZkUjraggCgpoLampO8c7yR988nsbh7D3mBRvJpYo5mMc3ppSvo4uiRFP4HK0lybiM9xcmcaemsnpNnYWHGPQsGM3JqKiAjH7JLPt12qRraGhJzYg23DnY3CKff34Gju5VyBV/AyRxnNZOkm7lymkvBDLyaTgfJfEtepMWfRr0/l45TbradSael8Ha+PPekBYQZtywYzMgZTTe51VRQ3djCsXde50yzjz/6F/FnSW+yzHEcl8NPreYCDpp0MqnSQwuTKJKLJEsvN0o9BTRxOljAb3pKmV/zFtvqHBYQZtyyYDAjprIrm99eXES5bzsZ2fPJrj/MFvdqHujKpjReRfRrJRyseJO3O27iHsd7fN/5HNNophsXTZpBmvi4qJNJo4cmzSRXWvFqGhc1k1R6yBEvHsd5vifPc1bzeL5nDQvOvsU7F1yU3JrGyraXQ60jCwgzDtioJDNiEn6TW7iVsG376xza+QJvtU/nHsd7LJCTzJB6VIRLmka+NJNCL256ecm/CjfdJNPDTDlHsZwnTXw0aA5BkkiSAB7Hef5X5/NM1SZ+1zGPwztf4A87Xqe6sWX8zQU1cM6rmgqo3hV5nDYH1rhjLQYzYuq8Pm5KaYq4yc2bks8n3pyR/+U1Few+0dDXSrjdsZ8/T/o9pXKKZPHjI4UuTSaDLl7130qRoxFBuc/5Hhc1k3Z146aHTOmgTnO5QZoI4MChDhwSBAlQJp9QSCOfqIcnzy9ixSsvQ2HZ2Ly8FG5Z1V7q5J3jjRxraCOlp5nMYCuLOU4DU2hnEoKykkO8xwJASKeTfJr5neMO0lLe+LSDPjvNWlBjmAWDGTEJu8ktHAp7du3gg54bWew4wSTtoizpKIqDszqVes3hcLCEqeKlx+Hi/2Mj5c5jNDk66UrKoEpvoMcfJODv5S7e55zmciQwneWO4xRKI6KKU3q50dFAqvaCA17tWcmd/S8vUTG6T4wDRmct9H1IXvAixbQyix5ypYVOUmnWdFY5DuAnCUHZ1lvGXa4PUcCJn6pgEQ/pb2npSMN3KplTp6BSUnhdVuFw2IiusUhUNdE1DEtZWZnu3bs30WWYQVTvf4NnDvcSyPSQkeqkrctPUksNm+e7KFmybmR+ab9QqOi5kcVygtl6lg1JH9JJKtXBaXyCh1cDK1nmrKJ20nz+wyI3y+bPG/yEFT55/vpoF5Nrd+IPKkcC0/ly0k4KHY241I9buunFydlgHsfx8Dq3scJ1kuW3rWdlPqPrm3O/MHjnHCzu2ks+TUyim5pALrc5D5EsfjLw0YYbP056cOLTZAqliTrNwYniR/qeT6IHt3T1vadD3ZzTHErkPF5No5tkFOgmhU+SZjDV1U3q1Jksmzsn1LK4bDT9PU0AIrJPVcsGfS1ewSAiTwL3ABdUdf4grwvwU+BuoBP4hqp+dK3PtWAY3ar3v8Hzn/j54/lJCMKS4ky+OU+Y5e6I7UmgpoIqXxov7a/DeewV9vpncbtjP7P1NOVyjHOaw/vBefwuuJI7kw6QJEJb8Rq+tHx6dCHVLyDSz75FUOFIYDr3Jr1PvjRRJBcQ4KxO5WBwJluDK1nr/JjpOZO46XN/QkleZuJOev0uE10OuLzgxSvCYBI91Gou6dLVd0IHSKUHt3RzIDiTxXKSAEKSaN/zTlLoJplUepgs7TTpZLLpGDQsbpAmkunFRYCLmkmPJIPAOUc+76XewS1T/SNzKery0Om2hr4bGuvP1XCo6gyHmx3k9J5Dg/E5F9YyjSxHB1McnTgdgj+oNAUn4Q2mUcx5HA7BGX74g0pvUK+orZZpeB1Zw2qNjZZgWA20A09fJRjuBr5DKBhuAX6qqrdc63MtGEa3qk8OUfnmM5zMW08g00NSSw0zG3dQunYzs25aEJtfEm4lHN3zFi92ldHtD/KQ/I7FHKNILrIvOIdKLWFncHF0rYRr/K59h4/y1IEOijuPsC84i02O3cyhhkJpJFn8+MOth2Ym82rgVmYkN4f7HlLje+09HAjVXZM49s4LeC9dJBAEX9DJ55IO4RoQBj51kSo9fSd2p/hxEdrnhN7AJLoJqIObHLV8ooUkoXRqKrMdtXSSQiepeDW9b0TXYGHhFD8ALj5tlXSqm2bNAKCJyaTRTbJTOO2cRUqwPeKkeb3amcRijnOSQgppJFNbmUw7FzSbQrlIL0kx+SuPRi8u6jSHQmnCRS+9uDinORSEQzPaz6jSQrocblKcSaSkZzNn43eG9P/UqAiGcCElwKtXCYZfAjtV9dnw80+ANapa/1mfacEwuv1kx3HEe4Zl7W/TkD6X/Paj7Eu/Hc2aznfXzxneh/cbirpn1w726yxWBveTrF3cKXuZRA+v6y2cJ489/lmUOU/iKC6LvpXwGar3v8ELe870XV7aGVjMvUm7mS8nKQ63Hk7qDaTQy/P+2ylIukR2ko/s7NyRb0EMCISzl3x83Ovhy0k7KXE04MLPWc0lo18YtGsqxXKRTlIivt2n0UWt5pIt7STTiyLs8s/jNudhBOjBhVfTmSpeqnVa6D3i+/SkP0hYXNY/NKbSQra0cYkMjgRLAJgtdQBDPmkOphcXzZrOjY4Gzmk2N8glzgWn4HFcoJMUenEiQDzOhgL04uSc5pBDK006mRukCZf4r/r7B9YmgB8XpyikyTmVXZmbmDt34ZD+n/qsYBhNnc+FwNl+z2vD264IBhF5GHgYwOPxxKU4c33qvD4KMj00MJei1o+onbyUQKaH+lhMv52RT9Wu5/nkxEUqehZyZ9JHzOckM+Us3bh4Qu8nVQIcccxhpbOK5bd9npWzY/NtvWTJOr6XU0FtTSq/PtrFsrN7eMW/EpJCJ8sbaMQjF6jXbL7jeolzmsOhwEyePD97ZEcv9buR7+wlH2/0hm7k+57rA1ClGydNmkG6+GjrFwYCfaOzJtENCDsCZbQ5Msh2dDKZdiB0CSPocvA8hX3f4L06iSp/PpuS3serGXg1g5QBI7pCw39baMNNb/i004uz756RZs3gkmYwWdrJkE7a1U0PTpp0MoVy8TNPmtFw0csUaaNesykQLx8Gb2KWnKNa8yPCKp5aSKciOJdSRw3J+IecSgL0ipPdLOKs5pEZwyntR1MwyCDbBv2rUtUtwBYItRhGsigzPIVZbsR7JmLIah15FGZNH94H11TwxzrhJ1VLuMe/jXsdu5ijZ5kuDdSTy2G9kbe0jFRnEl9M2cvc5etjFgp9POUUecr5nqeC3SdyyHv3JdoDafxj4EHudexmAScpkGa8mk6BNCNAbpKXI70lTKv5Pdvq+k2vMZxLTAP6ELR2L9W9eWxKep/vOyuZRjNBhIA4uaRpg4aBjxT8DhfPJt3DjNRO5kxL50+WzI+6ZVW9/w1278/jvYY2OnsCBIIKQT93ERkWudICRF5OuRwal8igUbMokoucpoDfB8uv+6Q5mMt9JXuCc1juOE6DZpGBj2xpG/6HD8HlFkOPOtnoqKBJJ5MrLUMKv8sthlYyWMnHBMVDYVZBzGocTcFQCxT3e14EnEtQLSZGNnm6aTjyDMfzNq38FDYAABMQSURBVHApcxl15DGzcQcLblwLNc3XfSLcfaKBQ7t20N6znHdZyPcdzzJbanlXF1FDPjuDS7jH9SGBm+5lzZKHYt/Z3Z+nnJVUkNJ9O08d6GBZ5xFeCa5EHRBUKHQ0Ua9TmCbNZEtb6OY6hH/vWUtBzVvsqut3iWkoQ1zDHeKvnQqQUv0mzp42pgFdQSePON/Dh4sZ0kA3LgI48auj776NYkcjnaQQcCSz1X0/qwqC3DMnj0euswVTsmTdlSFSU0Ftzdy++yJau/xM9nvxcJ5W3MxznKFJM2jSjL5LJTnShoMgM6mj2HGBep0y5JPmYC73ldRoLiscn3AuOIUbHKFpTrJoH9ZlqqHqGdDHkCx+6jSXAqK/XNbTr49hMp18MbCdOZ6bYlbjaOpj2Ah8m087nx9T1RXX+kzrYxjlwpc2Kvf8kV93LuNCUgFfzD7Nl1I/ZMr670dOsBfl5/UfjrpG9rNUjjFPTnNES/hIS9mpSyhPPhUaLhrrVsI19O97cAc7aGMSDYFsHnC+TVf4RN1BKhc0kxK5wDnNoSF8T0W6o4eD7hWsKghGjsyBTycjDLcO6s/VsPNYI+2NZ5gU7OB8MJPljuNMdzSQjJ/zmkm+XKIHF73q5KzmsTd4E9PES5tM4v2kpXwuX/nCgvy43ltQWd/CM++fofX4u6R2NRIIKsFgkIZAJpuSdjNNvBxjBov1GIXSSCepMeljcBGgVnOZKl5aNA0vGZwN5jHbUce+4Bwm04nDAQ7HYBcuYuuzRiUVBhsIEpqSwul0EAgogaDikMjaxtOopGeBNUAucB74EeACUNUnwsNVHwfuIjRc9Zuqes0zvgXD6FdZ38Jv39hNue8dLmXPJ/vSYSrcq3lg3cqhTcs9SCislEPMlDp+E1jDv+lGHkx6k2RnEjfd9sW4h0L/OmtrTvLaoQberRfmBY9TG8hhU9L7ZNNKvjTTKy661EmetFKjUxGUdnX3fYQA3ZJCpczAASzmOAeYQ6tjMnewh9xgM0eC0zml+axxfEy6+HCpn6A4+j63QadwTnPYGriV+Y4zIPQFwvryRSN3L8l1CF2KOswrjdNI67lIvlximqONXP+5mIxKqmUak+mglTSaCf2ba3ZMoWfSNG4v8LN63aZRM0V8Zf2AdUzmTxuR2kZFMIwUC4bR7yc7jtPi62Ve98d9HdBHUhaR6XZFP4piQCgskROs4AgrHJVsC95CFcW8ruWkpzj5P8q6r28o6ggYbPTSnyW9yVLHcZLxcyHc+VqnOWTSSba0XTHuX4BmTWemo4FTOo00umnUTIrkIiJwSdPJl2YAmjWtryVyUGeyNbCSpUlVHE2aMyoDwSTOWBmVZEar8M1jW2tS+r7FbPJ0R33dvs7rY3XvLvLbj3Em+9a+OZMaL3RH188wSCgs5ATljqO8FiynimL262zuT97DgvVfZ9mt17z9JW4Gjl5aV7uT9qCbd4OL8AWd3JJ0jDrNYbpc4Jzm9o3MuTzt9+Vx/5dH1JTIBV4OrCJH2gjgIFfbmCZNOAjSoy4y6OJnvV+kwHGJDOlkQ/JBWgpv5+9jMETXTBwWDOaaqnxpVL75DJK3noJMD+I9Q+WJHbB2M7OieH9pajPZjUfxaZCdF9x0di9i4/lfMScjGTL++rPffJVQuMOxn7eCS/pCoa9PoXAUtoD7jV6qrUn99BITx/mlfxObkt7HpX6KHedpIpMLmvXpuH+NvPv4reAiyh2VeDUNFDyOBnpwcjpYwN7gTUwVLzckedmdoD4EMz5YMJhr2lqTguStD92kRvgmtbz1VNWk8N0oBkKsLQrwz/UbOOftYj07OeGcQ08gyBGdydTPWpvhGqFwkNmRoZCoPoVohQPi4dtgw/43eGFPJotqd3IxmEUjWVQHbyCbVqZIW8S4f6Dv7uP5Uk2npjLXcYaTwRv4dXAdeXgJAOmOXirca1hVEOTvhzDU1JiBLBjMNQ33JrWSJetIPn2Qoo4KanzZ3Nh9gL2uJZxkKR1vv0/pyqwrT+jjLRQG6H+J6Z3jjRypa+FUZxr3JO3GoVwx7v+4FpIu3bSpm25S2B5cQY6jnW2OOxDgi9MusHFBPt+y1oGJAQsGc02xuEmtO6DMmpLC0rPv8K77DpZympRuJ+knK6la8nDkJalxHgp9wi2Ir9wWehoameNkR7+bxKbQggCTJTSixpuURYrTQa97GklZBXwjq5uFKzeMmhE1ZnywYDDXtMnTTeWJHewLT4R3+Sa10oWbo/6M0tRmkk7v56Pk5Xh6T3M4UMzd/t/wTNIXad9fx6PuDoC+GUwDZ/ewzz+LO5L2s1iOs0SqxlcoDGLQm8SMSQALBnNNs9wdsHYzVTUp1Ht9FGZNp3Th5tD2KK0tCvDXR8pJToL/2P3/cCcnORXIpyhwGj1Syb+eyybP0c5kbyUzA8k867+Te5N2c4scplgujvtQMGY0sfsYTNz87y8e5MjRQzzY+wKLtAoXvQRVeCewgBxpZbrjAmd0KgIUSiPFNJAh3WwLLOekeCwUjImhz7qPwRHvYszE9c15wuelgudcf8LTjnvp0BTc0sV9zveYl1TNGZ1KO25yxctMqcMtfn4bWE2VhYIxcWXBYOJmlruDjpkb6Q0oOcEmXtWVdOLGq+kEEDKkkxlSz2yppZtkPgzeTKe4OcgcCwVj4siCwcSPp5z7lxTyxZS9HHHMYaqjjbPBPCZJNx3qZrbUcrOc5bxO4Y3AMv4lcA8OgY2phywUjIkj63w2cTXL3cGF5XeQ9O5LuOmgiiLe7V3A153bmEwHnZrCab2BV4MrKXOepL34DjbOTaXIY6FgTLxYMJj46rd2weuHGqjpdrOsdx8ntZjfBabTmZzJIj3BfalHKLn1PmslGJMAFgwm/jzlLPOUM2lJCwd3b6e5IcBzPfdxOpCL2+UkdVoHDxY1WSvBmASxYDAJU1qQSemX/hcA/iLBtRhjPmWdz8YYYyJYMBhjjIlgwWCMMSaCBYMxxpgIFgzGGGMiWDAYY4yJYMFgjDEmQlyDQUTuEpFPRKRKRH4wyOtrRKRFRA6EHz+MZ33GGGPieIObiCQB/wysB2qBPSLyiqoeHbDru6p6T7zqMsYYEymeLYYVQJWqnlLVHuA54L44/n5jjDFRiGcwFAJn+z2vDW8b6FYR+VhEXheReYN9kIg8LCJ7RWRvY2PjSNRqjDETVjyDQQbZNnBd0Y+A6aq6CPgZ8PJgH6SqW1S1TFXL8vLyYlymMcZMbPEMhlqguN/zIuBc/x1UtVVV28M/vwa4RCQ3fiUaY4yJ5+yqe4DZIjIDqAMeBL7SfwcRyQfOq6qKyApCwdUUxxrNQDUVVPnS2FqTQp3XR2GWm02ebma5O2xKbGPGqbi1GFTVD3wb+ANQCfxGVY+IyCMi8kh4tweAwyLyMfAY8KCqDrzcZOKoypdG5ZvPIN4zFGSmIt4zVL75DFW+tESXZowZIXFdjyF8eei1Adue6Pfz48Dj8azJfLatNSlI3nqWtb9NA3PJbz/Kvrz1VNWk8N2bEl2dMWYk2J3P5jPVeX0EMj00pM+lqPUjGtLnEsj0UOf1Jbo0Y8wIsWAwn6kwy01SSw357UepnbyU/PajJLXUUJjlTnRpxpgRYkt7ms+0ydNN5Ykd7MtbH2opkMfMxh2ULtyc6NKMMSPEgsF8plnuDli7maqaFOq9PgqzplO6cHNouzFmXLJgMJ/NU84ssI5mYyYQ62MwxhgTwYLBGGNMBAsGY4wxESwYjDHGRLBgMMYYE8GCwRhjTAQLBmOMMREsGIwxxkSwG9zM1dlaDMZMSNZiMFdlazEYMzFZMJir2lqTwsnwWgyeln0sa3+bk3nr2VqTkujSjDEjyILBXJWtxWDMxGTBYK7K1mIwZmKyzmdzVbYWgzETkwWDuSpbi8GYickuJZnB2VBVYyYsazGYK9VUsPtEA0f3vMXepJX0ZhRT3PIRDUfegi88yqxE12eMGVFxDQYRuQv4KZAE/Iuq/njA6xJ+/W6gE/iGqn4U0yJqKqitOck7xxs51tBGZ0+AC8HJrHBUMsXRidMh+INKU3AS3mAaxZyP6a9PlFqmkeXoiOoYW4KpeDhPm07jq8nP0NWTyc3BE+wo+AuCNSm2mpsx41zcgkFEkoB/BtYDtcAeEXlFVY/22+0LwOzw4xbgF+E/Y6bKl8YnH7yNu+MSc/1BUoOdeOQCjf7JFEoTLnrpxcU5zaFAmkimN5a/PmF6cVHnz4nqGHvFRZ3mcLvjIOqHfK3kd6n38Z7/Jjw2VNWYcS+efQwrgCpVPaWqPcBzwH0D9rkPeFpDKoAsESmIZRFba1J42bmBVp1EBp3McdSSIZ0UOS5SzxQ6cFOvUyiUi6RILwjj4pEsvVEfo0t6KXRcJAhMEy9bg59De3twtZ21oarGTADxvJRUCJzt97yWK1sDg+1TCNT330lEHgYeBvB4PEMqos7ro1anslsWcxdvclEzAVCghXQqgnMpddSQjD+0cZyJ9hgzaadQmnghsBqfI5UjgZncEdjNGs/c+BZsjIm7eAaDDLJt4Gkpmn1Q1S3AFoCysrIhnb4Ls9y0nDvBSj1AMgFypQUXfnpw0qNONjoqaNLJoe3iHzfZIEBvlMfowk8mnRxVD04JckDnsEKqmLt8vQ1VNWYCiGcw1ALF/Z4XAeeuY59h2eTpZs7B7fRIJ21M4niwKNTHoJ/2MSSLnzrNpYDx08fQQ6jfILpjFPbqHJrIJsvhY2PyIUpuvY+Vs/NtqKoxE0A8g2EPMFtEZgB1wIPAVwbs8wrwbRF5jtBlphZVrSeGZrk7SL3l9itHJSVV0uS4IWLEzgfBReNrVFJSR1THWMs0vElZpKc4mTMtndVz8ijyWCgYM1HELRhU1S8i3wb+QGi46pOqekREHgm//gTwGqGhqlWEhqt+M+aFeMop8pTzldti/snGGDMuxPU+BlV9jdDJv/+2J/r9rMBfxbMmY4wxkWxKDGOMMREsGIwxxkSwYDDGGBPBgsEYY0wECfX3jl0i0gicGcJbcoGLI1TOaDYRj3siHjPYcU8kwznm6aqaN9gLYz4YhkpE9qpqWaLriLeJeNwT8ZjBjjvRdcTTSB2zXUoyxhgTwYLBGGNMhIkYDFsSXUCCTMTjnojHDHbcE8mIHPOE62Mwxhjz2SZii8EYY8xnsGAwxhgTYVwGg4jcJSKfiEiViPxgkNdFRB4Lv35QRJYmos5Yi+K4vxo+3oMisltEFiWizli71nH322+5iARE5IF41jdSojluEVkjIgdE5IiIvB3vGmMtin/jmSKyVUQ+Dh9z7GdojjMReVJELojI4au8HvvzmaqOqwehKb1PAjcCycDHwNwB+9wNvE5oYbNy4INE1x2n414JZId//sJEOe5++/2R0Oy+DyS67jj9984CjgKe8POpia47Dsf8t8B/C/+cBzQDyYmufZjHvRpYChy+yusxP5+NxxbDCqBKVU+pag/wHHDfgH3uA57WkAogS0QK4l1ojF3zuFV1t6peCj+tILRC3lgXzX9vgO8ALwAX4lncCIrmuL8CvKiqNQCqOtaPPZpjViBDRARIJxQM/viWGVuq+g6h47iamJ/PxmMwFAJn+z2vDW8b6j5jzVCP6T8S+pYx1l3zuEWkEPgi8ATjRzT/vecA2SKyU0T2icjX41bdyIjmmB8HSgktCXwI+E+qGoxPeQkT8/NZXBfqiRMZZNvAMbnR7DPWRH1MInIHoWAYD+vYRXPc/xP4a1UNhL5IjgvRHLcTWAasBdzA+yJSoarHR7q4ERLNMX8eOADcCcwEdojIu6raOtLFJVDMz2fjMRhqgeJ+z4sIfXsY6j5jTVTHJCILgX8BvqCqTXGqbSRFc9xlwHPhUMgF7hYRv6q+HJ8SR0S0/84vqmoH0CEi7wCLgLEaDNEc8zeBH2vo4nuViJwGbgY+jE+JCRHz89l4vJS0B5gtIjNEJBl4EHhlwD6vAF8P9+aXAy2qWh/vQmPsmsctIh7gRWDzGP7WONA1j1tVZ6hqiaqWAL8F/nKMhwJE9+/8d8DnRMQpIpOAW4DKONcZS9Eccw2hFhIiMg24CTgV1yrjL+bns3HXYlBVv4h8G/gDoVEMT6rqERF5JPz6E4RGptwNVAGdhL5ljGlRHvcPgRzg5+Fvz34d47NRRnnc4040x62qlSKyDTgIBIF/UdVBhzyOBVH+t/4vwFMicojQJZa/VtUxPRW3iDwLrAFyRaQW+BHggpE7n9mUGMYYYyKMx0tJxhhjhsGCwRhjTAQLBmOMMREsGIwxxkSwYDDGGBPBgsEYY0wECwZjjDERLBiMGWEikpToGowZinF357Mxo4GIPE9oxsslwJvAf01sRcZEz4LBmJGxAKhU1TsSXYgxQ2VTYhgTYyKSSmgytxtUdUwvEmMmJutjMCb25hFaXtFCwYxJFgzGxN4CQjOaGjMmWTAYE3sWDGZMsz4GY4wxEazFYIwxJoIFgzHGmAgWDMYYYyJYMBhjjIlgwWCMMSaCBYMxxpgIFgzGGGMi/P9n7qtvR6uNoQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(freud_rdf.bin_centers, freud_rdf.rdf, \"o\", label=\"freud\", alpha=0.5)\n",
    "ax.plot(*mdtraj_rdf, \"x\", label=\"mdtraj\", alpha=0.5)\n",
    "ax.set_xlabel(\"$r$\")\n",
    "ax.set_ylabel(\"$g(r)$\")\n",
    "ax.set_title(\"RDF\")\n",
    "ax.legend()"
   ]
  }
 ],
 "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.8.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}