results_report.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.signal import argrelextrema\n",
    "\n",
    "import simulate as sim\n",
    "import plotting as plot\n",
    "import observables as obs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Results for the final paper\n",
    "This file will contain things such as temperature or density loops so as to simulate and plot observables that give information about Argon's phase transitions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Temperature loop for the specific heat\n",
    "You may change the temperature range by changing the parameters:\n",
    "- temp_init\n",
    "- temp_final\n",
    "- temp_step\n",
    "\n",
    "And you should also set the other simulation parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input parameters\n",
    "particle_num = 4*(6**3) \n",
    "dim = 3 \n",
    "lattice_const = 1.5471 # [\\sigma]\n",
    "temperature = 5 # [\\epsilon/Kb]\n",
    "temperature_error = 0.01 # error in the temperature when rescaling [\\epsilon/KB]\n",
    "rescale_time = 0.1 # interval between rescalings\n",
    "\n",
    "run_time = 0.8 # sqrt(mass*\\sigma^2/\\epsilon)\n",
    "num_tsteps = 800 \n",
    "algorithm_method = \"verlet\" # options: \"verlet\" or \"euler\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "temp_init = 1.5\n",
    "temp_final = 2.5\n",
    "temp_step = 0.2\n",
    "\n",
    "temp_var = []\n",
    "c_variation = []\n",
    "c_error = []\n",
    "\n",
    "for temperature in np.linspace(temp_init, temp_final, int((temp_final-temp_init)/temp_step)+1):\n",
    "    print(\"Computing specific heat for temperature = \", temperature)\n",
    "    try:\n",
    "        init_vel = sim.init_velocity(particle_num, temperature)\n",
    "        init_pos, box_length = sim.fcc_lattice(particle_num, lattice_const)\n",
    "        eq_reached = sim.get_equilibrium(init_pos, init_vel, num_tsteps, run_time/num_tsteps, box_length, temperature, \"output_eq.csv\", method=algorithm_method, resc_thr=[temperature_error, rescale_time])\n",
    "        eq_pos, eq_vel = sim.load_final_data(\"output_eq.csv\")\n",
    "        temperature_eq = sim.temperature(eq_vel)\n",
    "        sim.simulate(eq_pos, eq_vel, num_tsteps, run_time/num_tsteps, box_length, \"output.csv\", method=algorithm_method)\n",
    "        c, Ac_autocorr, Ac_datablock = obs.specific_heat_error(\"output.csv\")\n",
    "        c_variation.append(c)\n",
    "        c_error.append(Ac_autocorr)\n",
    "        temp_var.append(temperature_eq)\n",
    "    except:\n",
    "        print(\"One step was skipped\")\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_xlabel(\"temperature $[\\epsilon/K_B]$\")\n",
    "ax.set_ylabel(\"$c_v$ $[K_B]$\")\n",
    "plt.errorbar(temp_var, c_variation, yerr = c_error ,marker = \".\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Serveral pair correlation functions\n",
    "Here we compute three pair correlation functions corresponding to solid liquid and gaseous phases. First you may set the simulation parameters and in each of the cases the required temperature:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Input parameters\n",
    "particle_num = 4*(5**3) \n",
    "dim = 3 \n",
    "lattice_const = 1.5471 # [\\sigma]\n",
    "rescale_time = 0.1 # interval between rescalings\n",
    "\n",
    "run_time = 0.8 # sqrt(mass*\\sigma^2/\\epsilon)\n",
    "num_tsteps = 500 \n",
    "algorithm_method = \"verlet\" # options: \"verlet\" or \"euler\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature = 0.1\n",
    "temperature_error = 0.01 \n",
    "init_pos, box_length = sim.fcc_lattice(particle_num, lattice_const)\n",
    "eq_reached = sim.get_equilibrium(init_pos, init_vel, num_tsteps, run_time/num_tsteps, box_length, temperature, \"output_eq.csv\", method=algorithm_method, resc_thr=[temperature_error, rescale_time])\n",
    "eq_pos, eq_vel = sim.load_final_data(\"output_eq.csv\")\n",
    "temperature_eq = sim.temperature(eq_vel)\n",
    "sim.simulate(eq_pos, eq_vel, num_tsteps, run_time/num_tsteps, box_length, \"output.csv\", method=algorithm_method)\n",
    "print(\"CALCULATING PAIR CORRELATION FUNCTION 1...\")\n",
    "r1, g1, eg1, eg11 = obs.pair_correlation_function_error(\"output.csv\", 0.01, box_length, r_max=3)\n",
    "positions = np.where(eg1>0.5)\n",
    "eg1[positions] = 0\n",
    "g1 = g1/g1.max()\n",
    "eg1 = eg1/g1.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature = 10\n",
    "temperature_error=0.5\n",
    "init_pos, box_length = sim.fcc_lattice(particle_num, lattice_const)\n",
    "eq_reached = sim.get_equilibrium(init_pos, init_vel, num_tsteps, run_time/num_tsteps, box_length, temperature, \"output_eq.csv\", method=algorithm_method, resc_thr=[temperature_error, rescale_time])\n",
    "eq_pos, eq_vel = sim.load_final_data(\"output_eq.csv\")\n",
    "temperature_eq = sim.temperature(eq_vel)\n",
    "sim.simulate(eq_pos, eq_vel, num_tsteps, run_time/num_tsteps, box_length, \"output.csv\", method=algorithm_method)\n",
    "print(\"CALCULATING PAIR CORRELATION FUNCTION 2...\")\n",
    "r2, g2, eg2,eg22 = obs.pair_correlation_function_error(\"output.csv\", 0.01, box_length, r_max=3)\n",
    "positions = np.where(eg2>0.5)\n",
    "eg2[positions] = 0\n",
    "g2 = g2/g2.max()\n",
    "eg2 = eg2/g2.max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature = 100\n",
    "temperature_error = 1\n",
    "init_pos, box_length = sim.fcc_lattice(particle_num, lattice_const)\n",
    "eq_reached = sim.get_equilibrium(init_pos, init_vel, num_tsteps, run_time/num_tsteps, box_length, temperature, \"output_eq.csv\", method=algorithm_method, resc_thr=[temperature_error, rescale_time])\n",
    "eq_pos, eq_vel = sim.load_final_data(\"output_eq.csv\")\n",
    "temperature_eq = sim.temperature(eq_vel)\n",
    "sim.simulate(eq_pos, eq_vel, num_tsteps, run_time/num_tsteps, box_length, \"output.csv\", method=algorithm_method)\n",
    "print(\"CALCULATING PAIR CORRELATION FUNCTION 3...\")\n",
    "r3, g3, eg3, eg33 = obs.pair_correlation_function_error(\"output.csv\", 0.01, box_length, r_max=3)\n",
    "positions = np.where(eg3>0.5)\n",
    "eg3[positions] = 0\n",
    "g3 = g3/g3.max()\n",
    "eg3 = eg3/g3.max()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we just plot them one after the other:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = [7.00, 6.00]\n",
    "plt.rcParams[\"figure.autolayout\"] = True\n",
    "fig, (ax1, ax2, ax3) = plt.subplots(3, sharex=True)\n",
    "fig.set_tight_layout(False)\n",
    "ax3.set_xlabel(\"$r/\\sigma$\")\n",
    "ax3.set_ylabel(\"$g(r)$\")\n",
    "ax2.set_ylabel(\"$g(r)$\")\n",
    "ax1.set_ylabel(\"$g(r)$\")\n",
    "\n",
    "ax1.errorbar(r1, g1, yerr=eg1, color = \"dodgerblue\", ecolor=\"lightgray\", label = \"T=0.1\")\n",
    "ax1.scatter(r1[positions],g1[positions], marker='o',color=\"white\")\n",
    "\n",
    "ax2.errorbar(r2, g2, yerr=eg2, color = \"goldenrod\",ecolor=\"lightgray\", label = \"T = 10\")\n",
    "\n",
    "ax3.errorbar(r3, g3, yerr=eg3, color = \"orangered\", ecolor=\"lightgray\", label=\"T=100\")\n",
    "\n",
    "ax1.legend()\n",
    "ax2.legend()\n",
    "ax3.legend()\n",
    "plt.savefig(\"results/W5_pair_correlation_functions_svg.svg\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Timing of the program including errors \n",
    "Here we will analyse a bit the efficiency of the code, now that we have observables and its errors. First, we analyse a loop where we compute the specific heat at different temperatures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def temperature_loop():\n",
    "    temp_init = 1\n",
    "    temp_final = 2\n",
    "    temp_step = 0.2\n",
    "\n",
    "    temp_var = []\n",
    "    c_variation = []\n",
    "    c_error = []\n",
    "\n",
    "    for temperature in np.linspace(temp_init, temp_final, int((temp_final-temp_init)/temp_step)+1):\n",
    "        print(\"Computing specific heat for temperature = \", temperature)\n",
    "        try:\n",
    "            init_vel = sim.init_velocity(particle_num, temperature)\n",
    "            init_pos, box_length = sim.fcc_lattice(particle_num, lattice_const)\n",
    "            eq_reached = sim.get_equilibrium(init_pos, init_vel, num_tsteps, run_time/num_tsteps, box_length, temperature, \"output_eq.csv\", method=algorithm_method, resc_thr=[temperature_error, rescale_time])\n",
    "            eq_pos, eq_vel = sim.load_final_data(\"output_eq.csv\")\n",
    "            temperature_eq = sim.temperature(eq_vel)\n",
    "            sim.simulate(eq_pos, eq_vel, num_tsteps, run_time/num_tsteps, box_length, \"output.csv\", method=algorithm_method)\n",
    "            c, Ac_autocorr, Ac_datablock = obs.specific_heat_error(\"output.csv\")\n",
    "            c_variation.append(c)\n",
    "            c_error.append(Ac_autocorr)\n",
    "            temp_var.append(temperature_eq)\n",
    "        except:\n",
    "            print(\"One step was skipped\")\n",
    "    return temp_var,c_variation,c_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext line_profiler"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%lprun -f temperature_loop temperature_loop()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusions\n",
    "\n",
    "As expected, in this function most of the time is spent in the functions sim.get_equilibrium and sim.simulate. This simulation in particular took 350 seconds, and the distibution looks like this:\n",
    "- simulate: 66.1%\n",
    "- get_equilibrium: 30.7%\n",
    "- specific_heat_error: 3.1%\n",
    "- rest: 0.1%\n",
    "\n",
    "The simulation involves 500 particles, a duration of 0.8 with 500 timesteps. We take them to equilibrium, time evolve them and compute the specific heat 5 times. Let's say we wanted to see the evolution of the specific heat with temperature in the three phases, we might want to make 10 measurements in each one, and that would take 35 minutes, or ~12 minutes per analysed phase. We regard this as plausible, though of course optimizing simulate would reduce the time cost of it. \n",
    "\n",
    "An option to reduce this is compute the specific heat while doing the simulation, instead of afterwards. We could do so by placing some booleans in the simulate, that would compute some parameters calling the specific heat function inside the simulate, where we already computed some of the quantities we know. However, as computing the specific heat only requires 3% of the time, we will leave it as it is.\n",
    "\n",
    "Another thing we could try though is to optimize the simulate function itself. That function spends most of the time computing the force acting on each particle. Thus, we tried reducing that cost by neglecting interaction between particles more than a certain distance away from each other. You still have to compute all the distances, then extract the indices of those that are far away... In the end it's not a lengthy operation due to its complexity, but due to the size of the operands, which we cannot reduce.\n",
    "\n",
    "__Now__ we do the same for the pair correlation function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_a_couple_pair_correlation():\n",
    "    try:\n",
    "        particle_num = 4*(5**3) \n",
    "        dim = 3 \n",
    "        lattice_const = 1.5471 # [\\sigma]\n",
    "        rescale_time = 0.1 # interval between rescalings\n",
    "\n",
    "        run_time = 0.8 # sqrt(mass*\\sigma^2/\\epsilon)\n",
    "        num_tsteps = 500 \n",
    "        algorithm_method = \"verlet\" # options: \"verlet\" or \"euler\"\n",
    "        temperature = 0.1\n",
    "        temperature_error = 0.01 \n",
    "        init_vel = sim.init_velocity(particle_num, temperature)\n",
    "        init_pos, box_length = sim.fcc_lattice(particle_num, lattice_const)\n",
    "        eq_reached = sim.get_equilibrium(init_pos, init_vel, num_tsteps, run_time/num_tsteps, box_length, temperature, \"output_eq.csv\", method=algorithm_method, resc_thr=[temperature_error, rescale_time])\n",
    "        eq_pos, eq_vel = sim.load_final_data(\"output_eq.csv\")\n",
    "        temperature_eq = sim.temperature(eq_vel)\n",
    "        sim.simulate(eq_pos, eq_vel, num_tsteps, run_time/num_tsteps, box_length, \"output.csv\", method=algorithm_method)\n",
    "        print(\"CALCULATING PAIR CORRELATION FUNCTION 1...\")\n",
    "        r1, g1, eg1, eg11 = obs.pair_correlation_function_error(\"output.csv\", 0.01, box_length, r_max=3)\n",
    "        positions = np.where(eg1>0.5)\n",
    "        eg1[positions] = 0\n",
    "        g1 = g1/g1.max()\n",
    "        eg1 = eg1/g1.max()\n",
    "    except:\n",
    "        pass\n",
    "    try:\n",
    "        temperature = 100\n",
    "        temperature_error = 1\n",
    "        init_vel = sim.init_velocity(particle_num, temperature)\n",
    "        init_pos, box_length = sim.fcc_lattice(particle_num, lattice_const)\n",
    "        eq_reached = sim.get_equilibrium(init_pos, init_vel, num_tsteps, run_time/num_tsteps, box_length, temperature, \"output_eq.csv\", method=algorithm_method, resc_thr=[temperature_error, rescale_time])\n",
    "        eq_pos, eq_vel = sim.load_final_data(\"output_eq.csv\")\n",
    "        temperature_eq = sim.temperature(eq_vel)\n",
    "        sim.simulate(eq_pos, eq_vel, num_tsteps, run_time/num_tsteps, box_length, \"output.csv\", method=algorithm_method)\n",
    "        print(\"CALCULATING PAIR CORRELATION FUNCTION 3...\")\n",
    "        r3, g3, eg3, eg33 = obs.pair_correlation_function_error(\"output.csv\", 0.01, box_length, r_max=3)\n",
    "        positions = np.where(eg3>0.5)\n",
    "        eg3[positions] = 0\n",
    "        g3 = g3/g3.max()\n",
    "        eg3 = eg3/g3.max()\n",
    "    except:\n",
    "        pass\n",
    "    return "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%lprun -f compute_a_couple_pair_correlation compute_a_couple_pair_correlation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusions\n",
    "In this case, it was the pair_correlation_function_error the one that took the most time. This simulation took 365 seconds and the distribution was approximately:\n",
    "- simulate: 20%\n",
    "- get_equilibrium: 20%\n",
    "- pair_correlation_function: 60%\n",
    "- rest: negligible\n",
    "\n",
    "So if we where to optimize anything in getting the pair correlation functions it would be the observable funcition itself. However, we are only running this code three times, on for each phase of Argon. So it's not worth our time to do so. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pressure loop at given temperature (Isotherms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0005951237015775439\n"
     ]
    }
   ],
   "source": [
    "# Input parameters\n",
    "particle_num = 4*(5**3) \n",
    "dim = 3 \n",
    "lattice_const = 1.5471 # [\\sigma]\n",
    "temperature = 2.74 # [\\epsilon/Kb]\n",
    "temperature_error = temperature/20 # error in the temperature when rescaling [\\epsilon/KB]\n",
    "rescale_time = 0.1 # interval between rescalings\n",
    "\n",
    "run_time = 1 # sqrt(mass*\\sigma^2/\\epsilon)\n",
    "num_tsteps = 500\n",
    "algorithm_method = \"verlet\" # options: \"verlet\" or \"euler\"\n",
    "\n",
    "MASS = 6.6335209E-26\n",
    "SIGMA = 3.405E-10\n",
    "RHO_conv = SIGMA**3/MASS\n",
    "print(RHO_conv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computing pressure for density =  0.5951237015775439  || lattice constant =  1.8871984946566494 || step:  1 / 11\n",
      "time step :     102/500 | temperature=2.63838/2.74000 (0.13023)\n",
      "time step :     500/500\n",
      "501/501\n",
      "Computing pressure for density =  0.6546360717352983  || lattice constant =  1.8281844882988239 || step:  2 / 11\n",
      "time step :     153/500 | temperature=2.73020/2.74000 (0.03027)\n",
      "time step :     500/500\n",
      "501/501\n",
      "Computing pressure for density =  0.7141484418930527  || lattice constant =  1.7759217771219482 || step:  3 / 11\n",
      "time step :     153/500 | temperature=2.68545/2.74000 (0.04885)\n",
      "time step :     500/500\n",
      "501/501\n",
      "Computing pressure for density =  0.7736608120508071  || lattice constant =  1.7291651099679997 || step:  4 / 11\n",
      "time step :     153/500 | temperature=2.68742/2.74000 (0.05891)\n",
      "time step :     500/500\n",
      "501/501\n",
      "Computing pressure for density =  0.8331731822085615  || lattice constant =  1.6869734058024999 || step:  5 / 11\n",
      "time step :     203/500 | temperature=2.66207/2.74000 (0.05468)\n",
      "time step :     153/500"
     ]
    }
   ],
   "source": [
    "density_init = RHO_conv*1000\n",
    "density_final = RHO_conv*2000\n",
    "density_step = RHO_conv*100\n",
    "\n",
    "density_var = []\n",
    "p_variation = []\n",
    "p_error = []\n",
    "beta_p = []\n",
    "beta_p_error = []\n",
    "b = int(((particle_num-1)/4)**(1/3)) + 1\n",
    "i = 1\n",
    "density_range = np.linspace(density_init, density_final, int((density_final-density_init)/density_step)+1)\n",
    "for density in density_range: \n",
    "    lattice_const = (particle_num/(b**3*density))**(1/3)\n",
    "    print(\"\\rComputing pressure for density = \", density,\" || lattice constant = \", lattice_const, \"|| step: \", i,\"/\",density_range.shape[0])\n",
    "    init_vel = sim.init_velocity(particle_num, temperature)\n",
    "    init_pos, box_length = sim.fcc_lattice(particle_num, lattice_const)\n",
    "    eq_reached = sim.get_equilibrium(init_pos, init_vel, num_tsteps, run_time/num_tsteps, box_length, temperature, \"output_eq.csv\", method=algorithm_method, resc_thr=[temperature_error, rescale_time])\n",
    "    eq_pos, eq_vel = sim.load_final_data(\"output_eq.csv\")\n",
    "    temperature_eq = sim.temperature(eq_vel)\n",
    "    sim.simulate(eq_pos, eq_vel, num_tsteps, run_time/num_tsteps, box_length, \"output.csv\", method=algorithm_method)\n",
    "    # compute pressure\n",
    "    P,P_autocorr,P_blocking = obs.pressure_error(\"output.csv\",temperature,box_length)\n",
    "    p_variation.append(P)\n",
    "    p_error.append(P_autocorr)\n",
    "    beta_p.append(P/temperature_eq)\n",
    "    beta_p_error.append(P_autocorr/temperature_eq)\n",
    "    density_var.append(density)\n",
    "    i += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Programs\\Anaconda3\\lib\\site-packages\\numpy\\core\\_asarray.py:83: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray\n",
      "  return array(a, dtype, copy=False, order=order)\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "setting an array element with a sequence.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;31mTypeError\u001b[0m: only size-1 arrays can be converted to Python scalars",
      "\nThe above exception was the direct cause of the following exception:\n",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-16-6d6ec7369aa6>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m      3\u001b[0m \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_xlabel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"density $[\\sigma^3/mass]$\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_ylabel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"$\\\\beta P$ $K_B/\\sigma^3$\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 5\u001b[1;33m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0merrorbar\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdensity_var\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbeta_p\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0myerr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mbeta_p_error\u001b[0m \u001b[1;33m,\u001b[0m\u001b[0mmarker\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m\".\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      6\u001b[0m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshow\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Programs\\Anaconda3\\lib\\site-packages\\matplotlib\\pyplot.py\u001b[0m in \u001b[0;36merrorbar\u001b[1;34m(x, y, yerr, xerr, fmt, ecolor, elinewidth, capsize, barsabove, lolims, uplims, xlolims, xuplims, errorevery, capthick, hold, data, **kwargs)\u001b[0m\n\u001b[0;32m   2984\u001b[0m                           \u001b[0mxlolims\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mxlolims\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mxuplims\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mxuplims\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2985\u001b[0m                           \u001b[0merrorevery\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0merrorevery\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcapthick\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcapthick\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 2986\u001b[1;33m                           **kwargs)\n\u001b[0m\u001b[0;32m   2987\u001b[0m     \u001b[1;32mfinally\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   2988\u001b[0m         \u001b[0max\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_hold\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mwashold\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Programs\\Anaconda3\\lib\\site-packages\\matplotlib\\__init__.py\u001b[0m in \u001b[0;36minner\u001b[1;34m(ax, *args, **kwargs)\u001b[0m\n\u001b[0;32m   1853\u001b[0m                         \u001b[1;34m\"the Matplotlib list!)\"\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mlabel_namer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1854\u001b[0m                         RuntimeWarning, stacklevel=2)\n\u001b[1;32m-> 1855\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0max\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1856\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1857\u001b[0m         inner.__doc__ = _add_data_doc(inner.__doc__,\n",
      "\u001b[1;32mC:\\Programs\\Anaconda3\\lib\\site-packages\\matplotlib\\axes\\_axes.py\u001b[0m in \u001b[0;36merrorbar\u001b[1;34m(self, x, y, yerr, xerr, fmt, ecolor, elinewidth, capsize, barsabove, lolims, uplims, xlolims, xuplims, errorevery, capthick, **kwargs)\u001b[0m\n\u001b[0;32m   3310\u001b[0m                 \u001b[0mxo\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0m_\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mxywhere\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlower\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnoylims\u001b[0m \u001b[1;33m&\u001b[0m \u001b[0meverymask\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   3311\u001b[0m                 \u001b[0mlo\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0muo\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mxywhere\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlower\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mupper\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnoylims\u001b[0m \u001b[1;33m&\u001b[0m \u001b[0meverymask\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 3312\u001b[1;33m                 \u001b[0mbarcols\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mvlines\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mxo\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlo\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0muo\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0meb_lines_style\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   3313\u001b[0m                 \u001b[1;32mif\u001b[0m \u001b[0mcapsize\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   3314\u001b[0m                     caplines.append(mlines.Line2D(xo, lo, marker='_',\n",
      "\u001b[1;32mC:\\Programs\\Anaconda3\\lib\\site-packages\\matplotlib\\__init__.py\u001b[0m in \u001b[0;36minner\u001b[1;34m(ax, *args, **kwargs)\u001b[0m\n\u001b[0;32m   1853\u001b[0m                         \u001b[1;34m\"the Matplotlib list!)\"\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mlabel_namer\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m__name__\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1854\u001b[0m                         RuntimeWarning, stacklevel=2)\n\u001b[1;32m-> 1855\u001b[1;33m             \u001b[1;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0max\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1856\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1857\u001b[0m         inner.__doc__ = _add_data_doc(inner.__doc__,\n",
      "\u001b[1;32mC:\\Programs\\Anaconda3\\lib\\site-packages\\matplotlib\\axes\\_axes.py\u001b[0m in \u001b[0;36mvlines\u001b[1;34m(self, x, ymin, ymax, colors, linestyles, label, **kwargs)\u001b[0m\n\u001b[0;32m   1051\u001b[0m                  for thisx, thisymin, thisymax in zip(x, ymin, ymax)]\n\u001b[0;32m   1052\u001b[0m         lines = mcoll.LineCollection(verts, colors=colors,\n\u001b[1;32m-> 1053\u001b[1;33m                                      linestyles=linestyles, label=label)\n\u001b[0m\u001b[0;32m   1054\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0madd_collection\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlines\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mautolim\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1055\u001b[0m         \u001b[0mlines\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mupdate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Programs\\Anaconda3\\lib\\site-packages\\matplotlib\\collections.py\u001b[0m in \u001b[0;36m__init__\u001b[1;34m(self, segments, linewidths, colors, antialiaseds, linestyles, offsets, transOffset, norm, cmap, pickradius, zorder, facecolors, **kwargs)\u001b[0m\n\u001b[0;32m   1252\u001b[0m             **kwargs)\n\u001b[0;32m   1253\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1254\u001b[1;33m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mset_segments\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msegments\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1255\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1256\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mset_segments\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msegments\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Programs\\Anaconda3\\lib\\site-packages\\matplotlib\\collections.py\u001b[0m in \u001b[0;36mset_segments\u001b[1;34m(self, segments)\u001b[0m\n\u001b[0;32m   1261\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mseg\u001b[0m \u001b[1;32min\u001b[0m \u001b[0msegments\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1262\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mseg\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mma\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mMaskedArray\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1263\u001b[1;33m                 \u001b[0mseg\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mseg\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfloat\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1264\u001b[0m             \u001b[0m_segments\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mseg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1265\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mC:\\Programs\\Anaconda3\\lib\\site-packages\\numpy\\core\\_asarray.py\u001b[0m in \u001b[0;36masarray\u001b[1;34m(a, dtype, order)\u001b[0m\n\u001b[0;32m     81\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     82\u001b[0m     \"\"\"\n\u001b[1;32m---> 83\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0morder\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0morder\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     84\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     85\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: setting an array element with a sequence."
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEQCAYAAABIqvhxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAFLtJREFUeJzt3X+w3XV95/HnS5DiDwqriTMuSQi7G4QUFfAWsLZbLbQN7GyYtmjJlAUsmoqiu0t1hh1bcHF2t+p23LrFYqgu2k5BZGpJaRSnLC6oRQm/f7MpomTilh8iXUX55Xv/+H6zOV7OzT2fJPece5PnY+bMfH98zve888m953W/vz7fVBWSJLV4waQLkCQtPIaHJKmZ4SFJamZ4SJKaGR6SpGaGhySp2djDI8mnkjyc5M4Z1ifJx5JsSnJ7kqPGXaMkafsmsedxCbBqO+tPAFb0r7XAn4yhJklSg7GHR1VdB3x3O01OAj5TnRuAA5K8cjzVSZJGMR/PeRwIPDQwv7lfJkmaJ/aedAFDZMiyoWOoJFlLd2iLl7zkJa879NBD57IuSdrt3HTTTY9W1eLW983H8NgMLB2YXwJsGdawqtYB6wCmpqZq48aNc1+dJO1GknxrR943Hw9brQdO66+6OhZ4oqq+M+miJEnbjH3PI8mlwBuBRUk2A+cDLwSoqouADcCJwCbgSeCt465RkrR9Yw+Pqlozy/oC3jWmciRJO2A+HraSJM1zhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWo2kfBIsirJfUk2JTl3yPplSa5NckuS25OcOIk6JUnDjT08kuwFXAicAKwE1iRZOa3Z7wGXV9WRwCnAx8dbpSRpeyax53E0sKmqHqiqp4HLgJOmtSngp/vp/YEtY6xPkjSLSYTHgcBDA/Ob+2WDPgCcmmQzsAF497ANJVmbZGOSjY888shc1CpJGmIS4ZEhy2ra/BrgkqpaApwI/FmS59VaVeuqaqqqphYvXjwHpUqShplEeGwGlg7ML+H5h6XOBC4HqKq/A/YFFo2lOknSrCYRHjcCK5IcnGQfuhPi66e1+TZwHECSw+jCw+NSkjRPjD08qupZ4GzgauAeuquq7kpyQZLVfbPfBd6e5DbgUuCMqpp+aEuSNCF7T+JDq2oD3YnwwWXnDUzfDbxh3HVJkkbjHeaSpGaGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR6SpGaGhySpmeEhSWpmeEiSmhkekqRmhockqZnhIUlqZnhIkpoZHpKkZoaHJKnZRMIjyaok9yXZlOTcGdq8JcndSe5K8hfjrlGSNLO9x/2BSfYCLgR+GdgM3JhkfVXdPdBmBfAfgDdU1eNJXjHuOiVJM5vEnsfRwKaqeqCqngYuA06a1ubtwIVV9ThAVT085holSdsxifA4EHhoYH5zv2zQIcAhSb6a5IYkq8ZWnSRpVmM/bAVkyLKaNr83sAJ4I7AEuD7J4VX1vZ/YULIWWAuwbNmyXV+pJGmoSex5bAaWDswvAbYMaXNlVT1TVd8E7qMLk59QVeuqaqqqphYvXjxnBUuSftIkwuNGYEWSg5PsA5wCrJ/W5q+ANwEkWUR3GOuBsVYpSZrRrOGR5B1JLk5ySpKrkpy1Mx9YVc8CZwNXA/cAl1fVXUkuSLK6b3Y18FiSu4FrgfdV1WM787mSpF1nlHMevwT8JnB9Vf18kot29kOragOwYdqy8wamCzinf0mS5plRDls91n+Zf6iff2oO65EkLQCjhMcfAVTVX/fzfzl35UiSFoJZw6Oq7t06neT8qvpfc1uSJGm+a73P4/wkLwZeBtwMXLb1LnBJ0p6j9VLdAn5EdzXUUuBrSV67y6uSJM1rrXse91bV+f30FUkuAS6iuyJLkrSHaN3zeDTJ67bOVNX9gLd2S9IepnXP4z3AZUluAu4AXgN8c5dXJUma10ba80jy+iSpqtuAI4BL+1XXAmvmqjhJ0vw06p7H6cCFSe4Hvgh8sar+Zu7KkiTNZyOFR1W9AyDJocAJwCVJ9qfb8/gi8NWqem7OqpQkzStNJ8yr6t6q+mhVraK7wuorwJuBr89FcZKk+WlnHga1HHimqt7d3zgoSdpD7MzzPP4EODzJeuDTSS7YRTVJkua5nQmPO6rqo8AjVfVmuiFLJEl7gJ0Jj2OT/DGwNMmrGf5scknSbmikcx5J9pp+NVVV/WySJcDr6E6aHzQH9UmS5qFRT5j/KfDW6QurajOwGbhyVxYlSZrfRj1stdewE+JJ9kvyuV1ckyRpnhs1PH4beH2St21d0A/FvhG4Zy4KkyTNX6PeYf5skl8HvpxkC3Ag8H7gbVX1t3NZoCRp/hn1hPnHgduBDwP/g25E3WOq6h/msDZJ0jw16gnzW+mGX3818ELgEODiJHfQ3e9x2RzVJ0mah0Y9bLVucL6/RHdrmPwrwPCQpD3IDo1tNXCJ7oZdW44kaSHYmTvMJUl7KMNDktRs1MfQnp7k0STfTfKZJPvNdWGSpPlr1D2P3wd+GTgU+Bbwn+esIknSvDfqCfN/rKpb+unfT+KTAyVpDzZqeLwyyVq6oUjupbvXQ5K0hxo1PM6nu6/jt+ju7Xhpkg3AbcDtVXXpHNUnSZqHdsVNgicChock7UFGvdrqbUmuSnJGkn2B0+geArWhqv7NnFYoSZp3Rr3a6r3AucCxwI10Y1v9A/Dfk5w+R7VJkuapUcPj6aq6E/h3wMHA7/SHsn4VeFfrhyZZleS+JJuSnLuddicnqSRTrZ8hSZo7o4bH55NcCZwAvLOqnuqXPwMsavnAJHsBF/bbWgmsSbJySLv9gPcAXhYsSfPMSOFRVefTfeEfB5yV5NtJrgGuAx5PcliSUYPoaGBTVT1QVU/Tjch70pB2H6R7fsiPRtyuJGlMRh5Vt6q+BHwJIEmAVwFHAkcAf9TPHzTCpg4EHhqY3wwcM9ggyZHA0qq6Ksl7Z9pQf+/JWoBly5aN+k+RJO2kHR2SvehuFryX9st0M2yT/39ltwfzUeCMEepYB6wDmJqaqlmaS5J2kUmMqrsZWDowvwTYMjC/H3A43fPSH6S7wmu9J80laf6YRHjcCKxIcnCSfYBTgPVbV1bVE1W1qKqWV9Vy4AZgdVVtnECtkqQhxh4eVfUscDZwNd1YWZdX1V1JLkiyetz1SJLa7dA5j51VVRuY9gjbqjpvhrZvHEdNkqTR+SRBSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUbCLhkWRVkvuSbEpy7pD15yS5O8ntSa5JctAk6pQkDTf28EiyF3AhcAKwEliTZOW0ZrcAU1X1GuAK4MPjrVKStD2T2PM4GthUVQ9U1dPAZcBJgw2q6tqqerKfvQFYMuYaJUnbMYnwOBB4aGB+c79sJmcCX5jTiiRJTfaewGdmyLIa2jA5FZgCfnGG9WuBtQDLli3bVfVJkmYxiT2PzcDSgfklwJbpjZIcD7wfWF1VTw3bUFWtq6qpqppavHjxnBQrSXq+SYTHjcCKJAcn2Qc4BVg/2CDJkcAn6ILj4QnUKEnajrGHR1U9C5wNXA3cA1xeVXcluSDJ6r7ZR4CXAp9LcmuS9TNsTpI0AZM450FVbQA2TFt23sD08WMvSpI0Mu8wlyQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTM8JAkNTM8JEnNDA9JUjPDQ5LUzPCQJDUzPCRJzQwPSVKziYRHklVJ7kuyKcm5Q9b/VJLP9uu/nmT5+KuUJM1k7OGRZC/gQuAEYCWwJsnKac3OBB6vqn8BfBT40HirlCRtzyT2PI4GNlXVA1X1NHAZcNK0NicBn+6nrwCOS5Ix1ihJ2o69J/CZBwIPDcxvBo6ZqU1VPZvkCeDlwKODjZKsBdb2s08luXNOKl54FjGtr/Zg9sU29sU29sU2r9qRN00iPIbtQdQOtKGq1gHrAJJsrKqpnS9v4bMvtrEvtrEvtrEvtkmycUfeN4nDVpuBpQPzS4AtM7VJsjewP/DdsVQnSZrVJMLjRmBFkoOT7AOcAqyf1mY9cHo/fTLwP6vqeXsekqTJGPthq/4cxtnA1cBewKeq6q4kFwAbq2o98Engz5JsotvjOGWETa+bs6IXHvtiG/tiG/tiG/timx3qi/gHvSSplXeYS5KaGR6SpGYLLjwc2mSbEfrinCR3J7k9yTVJDppEneMwW18MtDs5SSXZLS/THKUfkryl/7m4K8lfjLvGcRnh92NZkmuT3NL/jpw4iTrHIcmnkjw8071w6Xys76vbkxw160arasG86E6w/z3wz4B9gNuAldPavBO4qJ8+BfjspOueYF+8CXhxP33WntwXfbv9gOuAG4CpSdc9oZ+JFcAtwD/p518x6bon2BfrgLP66ZXAg5Ouew77418CRwF3zrD+ROALdPfYHQt8fbZtLrQ9D4c22WbWvqiqa6vqyX72Brp7anZHo/xcAHwQ+DDwo3EWN0aj9MPbgQur6nGAqnp4zDWOyyh9UcBP99P78/z7zXYbVXUd279X7iTgM9W5ATggySu3t82FFh7DhjY5cKY2VfUssHVok93NKH0x6Ey6vyx2R7P2RZIjgaVVddU4CxuzUX4mDgEOSfLVJDckWTW26sZrlL74AHBqks3ABuDd4yltXmr9PpnI8CQ7Y5cNbbIbGPnfmeRUYAr4xTmtaHK22xdJXkA3OvMZ4ypoQkb5mdib7tDVG+n2RK9PcnhVfW+Oaxu3UfpiDXBJVf1hktfT3Vt2eFX9eO7Lm3eavzcX2p6HQ5tsM0pfkOR44P3A6qp6aky1jdtsfbEfcDjw5SQP0h3TXb8bnjQf9ffjyqp6pqq+CdxHFya7m1H64kzgcoCq+jtgX7oBE/dEI32fDFpo4eHQJtvM2hf9oZpP0AXH7npsG2bpi6p6oqoWVdXyqlpOd/5ndVXt0IBw89govx9/RXchBUkW0R3GemCsVY7HKH3xbeA4gCSH0YXHI2Otcv5YD5zWX3V1LPBEVX1ne29YUIetau6GNllwRuyLjwAvBT7XXzPw7apaPbGi58iIfbHbG7EfrgZ+JcndwHPA+6rqsclVPTdG7IvfBS5O8u/pDtGcsZv+oUmSS+kOVS7qz/GcD7wQoKouojvncyKwCXgSeOus29xN+0qSNIcW2mErSdI8YHhIkpoZHpKkZoaHJKmZ4SFJamZ4SJKaGR7ShCQ5LMlFSa5Ictak65FaGB5a8JJ8IMl7d/E2v5bkgCTvbHzf8iQ/THLrbG2r6p6qegfwFrqxxwa384kkb2iresckeVGSW5M83d91Ls3K8JCGqKqfAw6gez5Mq7+vqiNGaZhkNfAV4Jppq46hG0ZlzlXVD/t6d9shybXrGR5akJK8v39K3N8CrxpYfmqSb/R/SX8iyV793sA9SS7un573pf6v7Zck+ZsktyW5M8lvDmzn+8AfAP+839ZHknwwyb8daPOfkrxnhFpfm+S6/ul9P073JMP/CFBV6/ug+q2B9ocB91fVc0k+l+SPk3wlybeS/HySzyS5P8knB95zcj/E+m1928X98tOT3NQ/He76mZZJzSb9hCtfvlpfwOuAO4AX0z3MZxPwXuAw4K+BF/btPg6cBiwHngWO6JdfDpwK/AZw8cB29x+Y/n7/vjsHli0Hbu6nX0D3pLqXT6tt+nv2Be4Fju7nP0g35ljoxhr6GN3gle8aeM85wG/30/cC5/TTF9CNgvtKuqfjPQ78VL/u5QPvPx94F91owncD+/TLDxi2bOB9DwKLJv3/62thvNzz0EL0C8Dnq+rJqvpHto2WehxdsNzYn3M4ju4xpADfrKqt5yFuovuSvwM4PsmHkvxCVT2xvQ+tqgeBx/rRin8FuKVmH1TweLrA+UY/fzvwsup8uareU1W/U1UXDrznV4EvJtmX7gv/v/XLfwh8sqq+U93T8Z4Enu7XndHvcd1Gd6jtR3QDH74I+MMkU9U9s2PYMqmZ4aGFatiIngE+XVVH9K9XVdUH+nWDzzJ5Dti7qu5n217Mf0ly3gif+6d0D5V6K/CpEdof3m9/q6OAm2dqnOTFdHsDW4CfoQuerQ8nei3w9b7dEmBLVVWS0+geu/pLVfVaur2Tu6p7BPHhwFeBdUneOWzZCP8G6XkMDy1E1wG/1p+32A/41/3ya4CTk7wCIMnLkhw000aS/FPgyar6c+C/0n2xD/q/dId5Bn0eWAX8LN1w37N5DHhN/3mHAL9O9zztmbwJuLaffjVw28C619DtuUAXJLcPtPtaVX0/yW8APwfckWRFVf2gqi4DrgL2HbZshH+D9DwL6nkeEkBV3Zzks8CtwLeA6/vldyf5PeBL6R49+wzdsf//M8OmXg18JMmP+7Y/ca9FVT2W7lnfdwJfqKr3VdXTSa4FvldVz41Q7qXA6n4bjwJrZjnUdQJwxUB93wDoD2G9qKoe79cNBsmngSuTnEz3XIYHquoH/UUFrwd+ANwFvB24aMgyqZnP85Aa9KF0M/DmqvrfQ9YvB66qqsN3cPs3A8dU1TM7U+cOfvaDwFRVPTruz9bC42EraURJVtJd2XXNsODoPQfsP8pNgsNU1VHjDo6tNwnSPVnux7O1l8A9D0nSDnDPQ5LUzPCQJDUzPCRJzQwPSVIzw0OS1MzwkCQ1MzwkSc0MD0lSM8NDktTs/wFihf60smJ02gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "ax.set_xlabel(\"density $[\\sigma^3/mass]$\")\n",
    "ax.set_ylabel(\"$\\\\beta P$ $K_B/\\sigma^3$\")\n",
    "plt.errorbar(density_var, beta_p, yerr = beta_p_error ,marker = \".\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "particle_num = 4*(5**3) \n",
    "density = 1000\n",
    "density_conv = RHO_conv*density\n",
    "lattice_constant = (particle_num/((((particle_num-1)/4)**(1/3)+1)**3)*density_conv)**(1/3)\n",
    "print(lattice_constant)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax.set_xlabel(\"density $[\\sigma^3/mass]$\")\n",
    "ax.set_ylabel(\"$P$ $\\epsilon/\\sigma^3$\")"
   ]
  }
 ],
 "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.6.6"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}