radial_distribution.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "1CzRnSUbYV2E"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "from tqdm import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "pi = math.pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "Q0fw2bTbYV2G"
   },
   "outputs": [],
   "source": [
    "rcut = 3.0  ##### Cut-off distance. If distance between atoms is greater than this, interaction between them is ignored.\n",
    "rcutsquare = rcut**2\n",
    "g2 = 1./rcutsquare\n",
    "g6 = g2*g2*g2\n",
    "potcut = 4.*g6*(g6-1) ##### Term to be subtracted from pair potential energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "i51VhivYYV2G"
   },
   "outputs": [],
   "source": [
    "##### Function to compute acceleration for a pair of atoms. r12 is relative displacement of atoms.\n",
    "\n",
    "def acceleration(r12): \n",
    "    r12square = np.dot(r12,r12)\n",
    "    f2 = 1./r12square #1/r^2\n",
    "    f6 = f2*f2*f2 #1/r^6\n",
    "    acc = 48.*f2*f6*(f6-0.5)*r12\n",
    "    return acc\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "id": "n8B7Y76QYV2H"
   },
   "outputs": [],
   "source": [
    "#### Function to compute potential energy for the system of particles. `pos' is the array of position vectors.\n",
    "def potentialenergy(pos):\n",
    "    potential = 0.\n",
    "    for i in range(Natoms-1):\n",
    "        for j in range(i+1,Natoms):\n",
    "            rij = pos[i] - pos[j] ## Relative position vector of the pair\n",
    "            for l in range(2):  ### Calculating the correct separation between particles because of periodic boundary conditions.\n",
    "                if abs(rij[l])>0.5*L: rij[l] -= L*np.sign(rij[l])\n",
    "            rijsquare = np.dot(rij,rij)\n",
    "            if rijsquare < rcutsquare: # Imposing interaction cut-off\n",
    "                f2 = 1./rijsquare\n",
    "                f6 = f2*f2*f2\n",
    "                potential += 4.*f6*(f6-1.) - potcut\n",
    "    return potential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "id": "tT49TSmdYV2H"
   },
   "outputs": [],
   "source": [
    "number = 5\n",
    "Natoms = number**2  # No. of atoms\n",
    "Ratom = 0.5  ### Radius of atom used to illustrate the atom as a sphere.\n",
    "rho = 0.25 ### Number Density\n",
    "L = math.pow(Natoms/rho,0.5) # Length of a side of the square containing the gas.\n",
    "T0 = 120. # Natural temperture scale, T0 = epsilon/k.\n",
    "# Temp = 0.01 #dimensionless temperature\n",
    "T = 2.5  # Temperature in Kelvin\n",
    "dt = 1e-2  # Time step for Verlet Algorithm\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {
    "id": "12feYni7YV2I"
   },
   "outputs": [],
   "source": [
    "###### Atoms are arranged initially in a square lattice configuration with spacing between neighbors equal to equilibrium separation of the potential energy.\n",
    "\n",
    "equilib_sep = math.pow(2.,1./6.)  ## Equilibrium separation of atoms r_0 = 2^{1/6}\n",
    "wall_spacing = (L-(number-1)*equilib_sep)/2.  ## Distance betwee the walls and the initial array of atoms.\n",
    "poslist = []  # List for positions of atoms\n",
    "vlist = []   # List for velocities\n",
    "\n",
    "################  Initialize positions and velocities ##############################################\n",
    "x, y = wall_spacing, wall_spacing\n",
    "\n",
    "for i in range(number):\n",
    "    for j in range(number):\n",
    "        xi, yi = x + equilib_sep*j, y + equilib_sep*i\n",
    "        poslist.append((xi,yi))\n",
    "\n",
    "\n",
    "####### Initially, all particles have the same speed but random velocity directions. ################\n",
    "v0 = math.sqrt(2.*T/T0)\n",
    "\n",
    "for i in range(Natoms):\n",
    "    phi = 2*np.pi*np.random.random()\n",
    "    vx = v0*np.cos(phi)\n",
    "    vy = v0*np.sin(phi)\n",
    "    vlist.append((vx,vy))\n",
    "###################################################################################################\n",
    "\n",
    "pos = np.array(poslist)  ### Converts lists to numpy arrays\n",
    "v = np.array(vlist) ## Scaled velocity\n",
    "\n",
    "\n",
    "V_cm = np.sum(v, axis = 0)/Natoms ## Correcting for CM velocity.\n",
    "V = np.array([V_cm,]*Natoms)\n",
    "v -= V  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta = 0.1 #bin size\n",
    "no_of_shells = int(0.5*L/delta)\n",
    "\n",
    "#RDF function\n",
    "def rdf(pair_distance, no_of_shells=no_of_shells):\n",
    "    r_array = np.zeros(no_of_shells)\n",
    "    which_bin = int(pair_distance/delta)\n",
    "    if which_bin<=no_of_shells-1:\n",
    "        bin_array[which_bin] += 1\n",
    "    for i in range(1,int(no_of_shells)):\n",
    "        r_array[i] = bin_array[i]/(2*rho*pi*i*delta*delta)\n",
    "    \n",
    "    return r_array  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "dWLHWFPEYV2J",
    "outputId": "5a0e698c-e84c-49fb-dd56-c1e70fd19277"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|███████████████████████████████████| 10000/10000 [00:51<00:00, 193.89it/s]\n"
     ]
    }
   ],
   "source": [
    "t_initial = 0. # Initial time.\n",
    "t_final = 100. # Time upto which simulation is carried out.\n",
    "rdf_data = []\n",
    "dist_array = np.zeros(no_of_shells)\n",
    "time = t_initial\n",
    "pair_distance_record = []\n",
    "counter = 0\n",
    "bin_array = np.zeros(no_of_shells)\n",
    "\n",
    "# potential_energy = potentialenergy(pos)\n",
    "kinetic_energy = 0.5*sum(np.square(v).sum(axis=1))\n",
    "# energy = kinetic_energy + potential_energy\n",
    "Time_List = [time]\n",
    "# Energy_List = [energy]\n",
    "# PotentialEnergy_List = [potential_energy]\n",
    "KineticEnergy_List = [kinetic_energy]\n",
    "# speed_list = [] #speed list for plotting histogram\n",
    "\n",
    "K_mean = T\n",
    "iterations = int((t_final-t_initial)/dt) # Number of iterations of the Verlet Algorithm.\n",
    "\n",
    "\n",
    "##################################### The Verlet Algorithm Loop ###################################\n",
    "\n",
    "for step in tqdm(range(iterations)):\n",
    "    pair_distance_list = []\n",
    "    pos += v*dt/2. ## Updating positions.\n",
    "    for i in range(Natoms): ## Imposing periodic boundary conditions.\n",
    "        for j in range(2):\n",
    "            if pos[i][j] > L: pos[i][j] -= L\n",
    "            elif pos[i][j] < 0: pos[i][j] += L\n",
    "    accel = np.zeros((Natoms,2)) ## Initialising acceleration array.\n",
    "    for i in range(Natoms-1):\n",
    "        for j in range(i+1,Natoms):\n",
    "            rij = pos[i] - pos[j]\n",
    "            for l in range(2):  ### Calculating the correct separation between particles.\n",
    "                if abs(rij[l])>0.5*L: rij[l] -= L*np.sign(rij[l])\n",
    "            if step>=2*iterations//5:\n",
    "                if step%10==0:\n",
    "                    counter += 1\n",
    "                    pair_distance = np.dot(rij, rij)\n",
    "                    dist_array += rdf(pair_distance)\n",
    "                    \n",
    "            if np.dot(rij,rij) < rcutsquare:\n",
    "                acc = acceleration(rij) # Computing acceleration for a pair of particles.\n",
    "                accel[i] += acc\n",
    "                accel[j] -= acc\n",
    "    v += accel*dt  ## Updating velocities.\n",
    "    pos += v*dt/2. ## Final updating of positions.\n",
    "    for i in range(Natoms): ## Imposing periodic boundary conditions.\n",
    "        for j in range(2):\n",
    "            if pos[i][j] > L: pos[i][j] -= L\n",
    "            elif pos[i][j] < 0: pos[i][j] += L\n",
    "\n",
    "    #         potential_energy = potentialenergy(pos)\n",
    "    kinetic_energy = 0.5*sum(np.square(v).sum(axis=1))\n",
    "\n",
    "    #scaling the velocity for keeping the temperature constant\n",
    "    if step>iterations//5:\n",
    "        Km = np.mean(kinetic_energy)/Natoms\n",
    "        lamda = np.sqrt(K_mean/Km)\n",
    "        v *= lamda\n",
    "    \n",
    "    #         energy = kinetic_energy + potential_energy\n",
    "    #         Energy_List.append(energy)\n",
    "    #         PotentialEnergy_List.append(potential_energy)\n",
    "#     KineticEnergy_List.append(kinetic_energy)\n",
    "#     time += dt\n",
    "#     Time_List.append(time)\n",
    "# #         if time>10:\n",
    "#             for atom in range(Natoms):\n",
    "#                 speeds = np.sqrt(np.dot(v[atom],v[atom]))\n",
    "#             speed_list.append(speeds)\n",
    "            \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "id": "BdOpCnOa3n0n"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x2049169cc70>]"
      ]
     },
     "execution_count": 83,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAk3ElEQVR4nO3deXzU1b3/8ddnJvsKIQtJAMMqBpDFiGuroli1KtbW7aqlvfrT29rFLlft7b0/b29t6+1me29rvVRtadXWXutCtVUrlbqgCCiogBhWWbKyZLJNkknO/SMTShFMSGbmOzN5Px8PHpn5JpnvZxDfOfl8z/ccc84hIiKJx+d1ASIiMjgKcBGRBKUAFxFJUApwEZEEpQAXEUlQKbE8WWFhoauoqIjlKUVEEt7q1asbnXNFhx6PaYBXVFSwatWqWJ5SRCThmdn2wx1XC0VEJEEpwEVEEpQCXEQkQSnARUQSlAJcRCRBKcBFRBKUAlxEJEEpwOPYe3vaePrtWrTkr4gczoAC3My2mdlbZrbGzFaFjxWY2Z/NrDr8cWR0Sx1+/vsv1fzTA6v50sNraO/s9rocEYkzRzMCP8s5N8s5VxV+fhuw1Dk3GVgafi4RtLupndz0FJ5Yu5uP/2w5O/a2eV2SiMSRobRQFgCLw48XA5cMuRr5O7VNQT40pZD7F57Ijn1tXPyTl1i+qdHrskQkTgw0wB3wrJmtNrMbwsdKnHM14ce1QMnhvtHMbjCzVWa2qqGhYYjlDi+1TUFK8jI4a2oxSz53OoU56Vxz3wrufXGL+uIiMuAAP905Nwc4H7jJzD588Cddb5ocNlGcc4ucc1XOuaqiovctpiVH0BzsorWzm9L8DADGF2bz2E2ncW7laO54agNf/t1aenoU4iLD2YAC3Dm3K/yxHngMmAvUmVkpQPhjfbSKHI5qm4IAjM7PPHAsJz2Fn10zh5vOmshjb+xi+eY9XpUnInGg3wA3s2wzy+17DJwLvA0sARaGv2wh8ES0ihyOavoCPC/j746bGZ+fN7n34uaaXV6UJiJxYiAj8BLgJTNbC7wGPOWcexq4E5hvZtXAOeHnEiG1gd4A72uhHCwj1c9Hpo/m6bdrCXZpeqHIcNXvhg7OuS3AzMMc3wOcHY2i5G8tlOK89MN+fsGsMh5ZvZNlG+s5b3ppLEsTkTihOzHjVG0gyKjsNNJT/If9/CkTRlGYk84Ta3bHuDIRiRcK8DjVN4XwSFL8Pi48vpSl79QTCHbFsDIRiRcK8DhV2xQ8bP/7YAtmldEZ6uGZt2tjVJWIxBMFeJyqDQQZ3U+Azxo7gnEFWSxZqzaKyHCkAI9Dwa5u9rZ2vm8K4aHMjAWzynh5UyP1zcEYVSci8UIBHofqAx0A/Y7AobeN0uPgqTdr+v1aEUkuCvA4VNPUDgwswCcV51JZmqfZKCLDkAI8Dn3QTTyHs2BWGWt27Gf7ntZoliUicUYBHocOtw7KB7loZhkASzQKFxlWFOBxqKYpSE56Cjnp/d4oC0DZiEzmji/g8TW7tMysyDCiAI9DdQOYQnioBbPK2NzQyvqaQJSqEpF4owCPQzVNwX6nEB7qgumlpPhMbRSRYUQBHocGMwIfmZ3GGVOKWLJ2tzZ6EBkmFOBxJtTdQ31zx1GPwAEWzC6npinIC9Xauk5kOFCAx5nGlk66e9xRj8ABzq0sYVxBFrcvWUd7p9YJF0l2CvA4c7RzwA+WkernPz9+PNv3tPH9ZzdGujQRiTMK8DhTG74L84OWkv0gp0wcxTUnj+P+l7eyevu+SJYmInFGAR5n+m7iGcwIvM9t5x9HWX4mtzyyVluuiSQxBXicqQkESfP7KMhOG/Rr5KSn8J1LZ7C5oZUfPVcdwepEJJ4owONMbVOQkvx0zGxIr/PhKUVcUTWWRS9sZu2O/ZEpTkTiigI8ztQ2BSnNG9gaKP35+oXHUZybwT8/spaOkFopIslGAR5nagNBSobQ/z5YXkYq3750Ou/WtfDTv2yKyGuKSPxQgMcR59yA9sI8GvOmlnDp7HLuXraZdbubIva6IuI9BXgc2d/WRUeoZ9BTCI/k/19USUaqnwdefS+irysi3lKAx5GaCEwhPJwRWWlML8/TSoUiSUYBHkfqwndhRnoEDlBZms87NQFC3T0Rf20R8YYCPI5EawQOMK0sj45QD1sbte2aSLJQgMeR2kAQn0FRbnrEX7uyLA9AbRSRJKIAjyO1Te0U5qST6o/8f5ZJxTmk+X2s260AF0kWCvA4UhvoiEr7BCDV72PK6BzWK8BFksaAA9zM/Gb2hpk9GX4+3sxWmNkmM3vYzAa/eIcAvSPwaFzA7DOtNJ91u5u08bFIkjiaEfgXgQ0HPf9P4C7n3CRgH3BdJAsbjmoifBPPoSrL8tjX1nVgzXERSWwDCnAzGwN8FLg3/NyAecAj4S9ZDFwShfqGjdaOEM3BEKPzI7MOyuFMC1/IXLdLbRSRZDDQEfiPgFuAvknEo4D9zrlQ+PlOoPxw32hmN5jZKjNb1dCgvRqPpG9UPDo/8jNQ+kwtzcNMM1FEkkW/AW5mFwL1zrnVgzmBc26Rc67KOVdVVFQ0mJcYFurCc8BHR2glwsPJSU+hYlS21kQRSRIpA/ia04CLzewCIAPIA34MjDCzlPAofAywK3plJr++m3gGs5nx0agszePNXfujeg4RiY1+R+DOua8558Y45yqAK4G/OOeuBp4HPhH+soXAE1Grchg40EKJ4iwU6L2QuWNvO03tXVE9j4hE31Dmgd8KfNnMNtHbE78vMiUNT7VNQfIzU8lM80f1PH13ZG5QH1wk4Q2khXKAc24ZsCz8eAswN/IlDU/RnkLY58BMlN0BTp4wKurnE5Ho0Z2YcaIuEIx6/xugODeDwpx03ZEpkgQU4HGipikY9f53n2lleZqJIpIEFOBxoDPUQ2NLR0xG4NDbB99U36KNjkUSnAI8DtQ3x2YGSp9pZXmEehzVdS0xOZ+IRIcCPA7UxmgOeJ/K0vDa4OqDiyQ0BXgc6JsDXhrFdVAOVjEqm6w0v26pF0lwCvA4cGAEHqMWis9nHFeqC5kiiU4BHgdqmoJkpvrJyzyqaflDMq0sjw01zfT0aG1wkUSlAI8DteE54L2r9MZGZWkeLR0h3tvbFrNzikhkKcDjQG1TkJK86C0jezjTyvIBLS0rksgU4HFgX2sno3JiG+CTS3Lw+0x9cJEEpgCPA4FgF/mZqTE9Z0aqn8nF2uRYJJEpwD3mnCPQHiIvI7YBDr198HUKcJGEpQD3WEeoh87unpjOQOlTWZZHfXMHDc0dMT+3iAydAtxjfRsreDICDy8tqwuZIolJAe6xQF+Ax7gHDjCtNDwTRW0UkYSkAPdYINgb4LG+iAmQn5VK+YhM3ty5P+bnFpGhU4B7LNAeAiAvI/Y9cID5lSU8va6WZRvrPTm/iAyeAtxjfSNwL1ooALeeN5VjS3K5+eE17NBdmSIJRQHuMS8vYgJkpvm555oT6O52fPbB1wl2aZMHkUShAPfY3y5ietNCAagozOYHl8/krV1N/MeT6z2rQ0SOjgLcY4FgiIxUH+kpfk/rOHfaaP7pjIk8tOI9Hlm909NaRGRgFOAeC7R3edY+OdRXz53CyRMK+Ppjb2lqoUgCUIB7rKm9y7MLmIdK8fv476vmkJ+ZymceXH2gPy8i8UkB7jEvFrL6IEW56dx99Rx27Wvnyw+voau7x+uSROQIFOAe613IyrsLmIdTVVHA7RdVsvSdej730Ot0hhTiIvFIAe6xQDB+WigHu/aUCm6/qJJn1tXxTw+s1vRCkTikAPdYUxxdxDzUp08bz7c/NoPnN9Zz/eJVtHWGvC5JRA6iAPdQ71rg8dUDP9Q/nDSO739iJss3N/Kp+1fS0qEQF4kXCnAPtXZ20+O8vYlnID5+whh+fOVsVr+3j2vvW6HZKSJxot8AN7MMM3vNzNaa2Toz+0b4+HgzW2Fmm8zsYTNLi365ySXg8W30R+OimWXcffUc3t7VxNX3vqoQF4kDAxmBdwDznHMzgVnAeWZ2MvCfwF3OuUnAPuC6qFWZpJo8XAt8MD4ybTSLPlnFxtpmbvjVKjpCurAp4qV+A9z1agk/TQ3/ccA84JHw8cXAJdEoMJkl0gi8z1nHFvP9y2ayYutevvK7tfT0OK9LEhm2BtQDNzO/ma0B6oE/A5uB/c65vitaO4HyI3zvDWa2ysxWNTQ0RKDk5BEI9v71xfNFzMNZMKucr50/lSffrOE7f9rgdTkiw9aAAtw51+2cmwWMAeYCUwd6AufcIudclXOuqqioaHBVJql4WIlwsG748AQ+dWoFP39xK/e/tNXrckSGpaNKDufcfjN7HjgFGGFmKeFR+BhgVzQKTGYHNnNIoBZKHzPj3y6spKapnW8+tZ7R+RlcMKPU67JEhpWBzEIpMrMR4ceZwHxgA/A88Inwly0EnohSjUmr7yJmbpzdSj9Qfp/x4ytnc8K4kdz88Bpe27rX65JEhpWBtFBKgefN7E1gJfBn59yTwK3Al81sEzAKuC96ZSanQHuInPQUUvyJOx0/I9XPzz9ZxZiRmVy/eCXbGlu9Lklk2BjILJQ3nXOznXPHO+emO+f+I3x8i3NurnNuknPuMudcR/TLTS6BYFfcLWQ1GCOz01j86bkEQz38cvk2r8sRGTYSd+iXBAJxtBb4UI0tyGLescU8+WYN3ZpaKBITCnAPxfNCVoNx8awyGls6eGXzHq9LERkWFOAeCgRDSTMCB5g3tZic9BSWrNWEJJFYUIB7qLeFkvg98D4ZqX7OnVbCn96u1W32IjGgAPdQ70XM5BmBA1w8s4zmYIhlG3XXrUi0KcA90t3jaE6yFgrAaZMKGZWdxpK1u70uRSTpKcA90pKg66D0J9Xv44IZpSzdUKfNH0SiTAHukb/dRp88PfA+F88qI9jVw5/X13pdikhSU4B7JNHWAj8aJ4wbSVl+BkvWqI0iEk0KcI8k4lrgA+XzGRfNKuPF6kb2tnZ6XY5I0lKAe+RACyWJphEe7OKZZYR6HH98q8brUkSSlgLcI4H25LyI2aeyNI+JRdmajSISRQpwj/xtBJ6cAW5mLJhVzspte9m9v93rckSSkgLcI4H2LswgJy05WyjQ20ZxDp58U6NwkWhQgHukqb2L3PQUfD7zupSoqSjM5vgx+WqjiESJAtwjgWCI/KzkbJ8c7OKZZby9K8DmhhavSxFJOgpwjwSSbCnZI7loZhlmsFgbPYhEnALcI8m4kNXhlORlcNXccfzqle389PlNXpcjklSS9wpanGtq72J8YbbXZcTENxdMp60jxPee2UiKz7jxjIlelySSFBTgHgm0h5J2Dvih/D7j+5fNJNTj+M6f3sHvM67/0ASvyxJJeApwjwyXFkqfFL+Pu66YRXeP446nNpDq97Hw1AqvyxJJaOqBe6Cru4e2zu6kvYnnSFL9Pv7rqtnMryzh9iXreODV7V6XJJLQFOAe+NtCVsPvF6BUv4+f/MNs5k0t5l8ff5vfrdzhdUkiCUsB7oFA32YOw2Ae+OGkp/i5++o5fGhyIf/y2Fu8tnWv1yWJJCQFuAeSeSnZgcpI9fPTq+cwtiCLzz74OnWBoNcliSQcBbgHkn0hq4HKy0jlnmtOoLUjxGcffJ3OUI/XJYkkFAW4B5o0Aj/g2NG5fPcTx7N6+z6+9dR6r8sRSSgKcA/0rQWerJs5HK2LZpZx3enjWfzKdh57Y6fX5YgkDAW4B/paKMPlRp6BuO38qZw0voCvPfoW63cHvC5HJCEowD0QaO8ixWdkpvq9LiVu9E4vnMOIzDRufGAV+9u0l6ZIf/oNcDMba2bPm9l6M1tnZl8MHy8wsz+bWXX448jol5scAsEu8jJTMUvetcAHoyg3nbuvmUNtU5CbHnqdnfvavC5JJK4NZAQeAr7inKsETgZuMrNK4DZgqXNuMrA0/FwGoKk9NCxv4hmIOeNG8q1LZvDa1r2c+b1lfPV/12otcZEj6DdFnHM1QE34cbOZbQDKgQXAmeEvWwwsA26NSpVJJtDepf73B7j8xLGcPrmQRS9s4bcr3+P3r+/kgumlfPasiUwry/e6PJG4cVQ9cDOrAGYDK4CScLgD1AIlR/ieG8xslZmtamhoGEqtSaOvhSJHVjYik3+/eBov3TqPz5wxkRfebeCj//US1y9eeWAapshwN+AAN7Mc4PfAzc65v5sm4JxzgDvc9znnFjnnqpxzVUVFRUMqNlkMl914IqEwJ51bzpvKS7fN4yvzp/CXd+r5wbMbvS5LJC4MKMDNLJXe8H7QOfdo+HCdmZWGP18K1EenxOTT1B7SHPCjlJ+ZyufPnswnT6nggVe3a6qhCAObhWLAfcAG59wPD/rUEmBh+PFC4InIl5ec1EIZvC+dM4URWWn8+5J19P7iJzJ8DWQEfhpwLTDPzNaE/1wA3AnMN7Nq4Jzwc+lHsKubzlCPWiiDlJ+Vyq3nHctr2/ayZO1ur8sR8dRAZqG8BBxpwvLZkS0n+Wkhq6G77ISxPLTiPb711AbOPq6EnHS1o2R40p2YMTacN3OIFJ/P+MaC6dQ3d/DfS6u9LkfEMwrwGGs6sJCVRuBDMWvsCC6vGsP9L29lU71u9JHhSQEeY1rIKnJuOW8qGal+vvEHXdCU4UkBHmPajSdyCnPS+fL8KbxY3cgz6+q8Lkck5tSIjbG+/TA1Dzwyrj35GH772g6++eR6CnPSaGzppKE5SENzB/XNHTS2dHDBjFIunTPG61JFIk4pEmMagUdWit/HNxZM48pFr/KJe145cNxnMConnRSfsfSdekZkpTJv6mFXexBJWArwGAu0d5Ge4iNDa4FHzMkTRvHQ9SfREeqhKDed4tx0CrLTSPH7aO/s5rL/Wc4XfrOGRz97KlNKcr0uVyRi1AOPMd2FGR2nTirkrKnFTC/PpzgvgxR/7z/tzDQ/P/9kFZlpfq5fvIq9rdooQpKHAjzGAloLPOZK8zNZdO0J1AaCfOaB1XSGerwuSSQiFOAx1tSuEbgXZo8byXc/fjwrtu7ldq2jIklCQ8EYCwS7KMhO87qMYemS2eW8W9fM3cs2M3V0LgtPrfC6JJEh0Qg8xrQWuLe+eu6xzK8s4T+eXM9f39UGI5LYFOAxFghqLXAv+XzGXVfMYnJxDgvvf41r71vBc+vr6O5RS0USjwI8hpxzvT1wjcA9lZOewm9vOJmvzJ9CdV0L1/9qFWd873nu+etm9mmWiiQQBXgMtXV2093jtA5KHBiRlcbnz57Mi7eexc+unsOYkZnc+ad3OOk7S/nn/13Lzn1tXpco0i/9Lh9DWgs8/qT6fZw/o5TzZ5SysbaZX7+6jf9dtZMn1u7mutPH89kzJ5Kr35gkTmkEHkOBvqVkFQhx6djRudxxyQye/+qZfHRGKT9btpkzv7eMB17dTqhbc8cl/mgEHkNNfeug6CJmXCsbkcldV8zi06dVcMdTG/jXx99m8fJtfPUjxzJ2ZBad3T10hnr/dHX3EOpxnDShQD+YJeaUJDGkhawSy/FjRvDwDSfz7Po67vzTO9z469VH/NpxBVncu7BKa61ITCnAY0ibOSQeM+Mj00Yzb2oxL7zbQFd3D2kpPtL8ftJSfKT6jb2tndz26Ft87Kcv88MrZvGRaaO9LluGCQV4DB0YgSvAE06q38fZxx15Odo/fC6fG3+9iht/vZovnTOFz8+bhM93pL3ARSJDFzFjqG8zh1wtZpV0Rudn8PCNp3Dp7HLueu5dbnrodVo7Ql6XJUlOSRJDTe1dZKX5SfXr52Yyykj184PLZ1JZlse3/7iBrY2tfOfSGRxXmqf13yUqFOAxFGjvUv87yZkZ139oAlNKcvncQ6/zsbuX47Pei5yTS3KZXJzDlJJcppTkMqk4h7QU/TCXwVOAx1AgqNvoh4sPTyniua+cwWtb91Jd10J1fTPVdS08/049ofC6Kyk+Y1JxDseV5jF1dC7HleYxZmQmNU1Btu9pY/ueVrbtaWX7njZqmoJUjMri+DEjOH5MPjPHjmBiUQ5+9dmHNQV4DAXatZDVcFKcm8GFx5f93bHOUA/b9rSysbaZDTUBNtQEeGXzHh57Y9f7vj/N72NsQSYVo7KZc8xItja08tgbu/j1q9sByE7zM708nxMrCjh14ijmHDNSrZphRmkSQ03tXZTmZ3hdhngoLcV3oIVy0cy/hfu+1k421ATYtb+d8hGZHFOYzei8jPeNsHt6HFsaW1i7o4k3d+5nzc4mfvbXzfzk+U2kpfioOmYkp00q5JSJozi+PP/A1nKSnBTgMRQIdjF1tG70kPcbmZ3GqZMK+/06n8+YVJzLpOJcPn7CGACag12s3LaXlzftYfnmPXzvmY0AjBmZyS3nTeWi40sxU6slGSnAYyig7dQkCnIzUpk3tYR5U3vnqe9t7eSlTY3cs2wzX/jNG9z30la+fsFxzB1f8IGv45xT0CcYBXiMdPc4mjtCmgMuUVeQncbFM8u4cEYpj76xi+8/s5HL/+cVzq0s4bbzpzKhKAfnHFsaW1m9bR+rtu9l9fZ97NjbzkePL+X6D41nWlm+129DBqDfNDGz+4ELgXrn3PTwsQLgYaAC2AZc7pzbF70yE9+ufe041/trrUgs+HzGJ04Yw0dnlHLfS1v42bLNnHvXC5xYUcA7tQH2tf1taYcTjhlJ1TEF/OHN3Tz2xi5OnTiK//ehCZwxpUh3lMYx6293bjP7MNAC/OqgAP8usNc5d6eZ3QaMdM7d2t/Jqqqq3KpVqyJQduJZuqGO6xav4vefOZUTjhnpdTkyDDU0d/Cj595l9fZ9zCjP7w3tipFMKMw5ENJN7V385rX3+OXL26gNBJlUnMP1p4/nktnlmuHiITNb7Zyret/x/gI8/M0VwJMHBfhG4EznXI2ZlQLLnHPH9vc6wznA7/nrZu780zusvf1c3cwjca8z1MNTb+3m5y9sZX1NgJFZqVw1dxzXnHwMZSP0W2SsHSnAB9uQLXHO1YQf1wJHXOXHzG4AbgAYN27cIE+X+KrrWijJS1d4S0JIS/HxsdljuGRWOa9s2cPi5du456+b+Z8XtvCRaSV86tTxnFgxUhc9PTbkK2rOOWdmRxzGO+cWAYugdwQ+1PMlqk31zUwqzvG6DJGjYmacOrGQUycWsmNvGw+8up3frtzBH9+qpbI0j+nlebR2dtPe2U1bZyj8sZtjRmVzy3nHDmh99N3723l8zS6c611PJiPVR2aqn4xUP1lpfuaOLyArTRf/D2ewfyt1ZlZ6UAulPpJFJRvnHNX1LVxeNdbrUkQGbWxBFl+74DhuPmcKj6/ZxQOvbueFdxvJSvOTmdYbtvlZaYzO9/HK5j2c/+MX+Ye54/jS/CkUZKe97/Xqm4Pc/fxmHlrxHp0fsGVdWX4Gt188jXMrSzTiP8RgA3wJsBC4M/zxiYhVlIR2NwVp6+zWCFySQmaan6vmjuOquUduie5r7eSu597lwRXv8cSaXdx8zhSuPeUYUv0+9rV2cs8Lm1m8fBtd3Y7Lq8Zw01mTKMxJp6Orh/auboJd3QRD3eze3853n97Ijb9ezdlTi/n3i6cxtiDrfefrDPXw5/V1/Hble2yoCRwYvWem9v5wyUz1U5ybwcJTK6gsy4vmX09MDWQWym+AM4FCoA64HXgc+B0wDthO7zTCvf2dbLhexFy2sZ5P/WIlD99wMidNGOV1OSIx825dM998cj0vVjcyoSibeccW89uVO2jtDHHJrHK+ePZkKgqzP/A1urp7+OXL27jruXfpcY4vnj2F604fT1qKjy0NLTy8cgePrN7JntZOykdkctqkUYS6HW2d3bR39bZ32ru62dbYSnNHiPOmjebm+ZOZOjpxgnxIs1AiZbgG+L0vbuGOpzbw+r/NP+yvkiLJzDnH8xvruePJDWxpbOX86aP50vwpR71/6O797XzjD+t4Zl0dk4tzKMxJ55Ute/D7jHOOK+bKueP48OSiI67Q2NTexX0vbeUXL22luSPER2eU8sVzJifEPqYKcA/d+sibPLehjtX/Nt/rUkQ809Xdw/62Lopy04f0Oks31HHHUxvo7nFcceJYLjthDMV5A18krqmti3tf2sIvXt5Ga2dvkH/ylIqozappbOnggVe384V5kwd9U1SkpxHKUajWDBQRUv2+IYc3wNnHlXzg/qT9yc9K5SvnHss/njaee1/awuLl23nyzRrGFWRx6ZxyPj5nzPv67J2hHtbs2M/Lmxp5ZcseinLT+dI5k5lUfOTRu3OOR1bv5Ft/3EBrR4izji1m5tgRg677cDQCjzLnHDO/8SwXzyrjjktmeF2OiByirTPE02/X8vvXd7J88x6cg7njC/jY7HJagiFe3tzIa1v30tbZjc9gWlk+WxtbaesMccWJY7n5nCmUHPIbwLbGVv7lsbdYvnkPVceM5M6Pz/jAsO+PRuAeaWjuIBAMMXkI//FEJHqy0lK4dM4YLp0zhl3723n8jV38fvVOvvboWwBMLMrmEyeM4bRJhZw8YRT5mansbe3kJ3/ZxK9f3cZjb+ziutPHc+MZE8lM9bPohS3819Jq0vw+vvWx6Vx14riorSejEXiUvbypkavvXcFD1580oPWeRcR7zjneqW1mZFYaoz9gE5Yde9v4wbMbeXzNbkZkpVKYk86m+hYumDGa2y+a9r6R+WBpBO6R6rpmAPXARRKImXFcaf/TDMcWZPGjK2dz/Ycm8N1nNrJ7fzs//2QV8ysH36M/GgrwKKuubyEvIyUiF29EJD5NL8/nV/84N+bn1YZ5UVZd38LkklzdAiwiEacAj7JN9S1MVvtERKJAAR5Fe1o62Nvaqf63iESFAjyKqutbAJicALfqikjiUYBH0YEA1whcRKJAAR5Fm+tbyE7zU/oB80hFRAZLAR5F1fXNTNIMFBGJEgV4FFXXaQaKiESPAjxKmtq6qG/u0AwUEYkaBXiUbGrovYVeI3ARiRYFeJRU1/XNQNEUQhGJDgV4lFTXt5CR6qN8ZKbXpYhIklKAR0l1fQsTi3KOuD+fiMhQKcCjZFNds/rfIhJVCvAoaA52sbspqFvoRSSqFOBRsLmhFdAmDiISXQrwKNikNVBEJAYU4FFQXd9Mmt/HuIIsr0sRkSSmAI+CTXUtjC/MJsWvv14RiR4lTBRU17cwqUTtExGJLgV4hLV3drNjX5v63yISdQrwCOrpcdz/8lac0y30IhJ9KV4XkCy2NrZyyyNrWbltH/OmFnP2ccVelyQiSW5IAW5m5wE/BvzAvc65OyNSVQLp7nH8cvk2vvfMO6T5ffzgsplcOqdcmziISNQNOsDNzA/8FJgP7ARWmtkS59z6SBUX77Y0tHDLI2+yavs+zp5azLcvnUFJnrZPE5HYGMoIfC6wyTm3BcDMfgssACIe4F9/7C1e27o30i87ZO/tbSM9xccPL5/Jx2Zr1C0isTWUAC8Hdhz0fCdw0qFfZGY3ADcAjBs3blAnKhuRyeQ4nJZ30oQCvjBvMsUadYuIB6J+EdM5twhYBFBVVeUG8xo3nTUpojWJiCSDoUwj3AWMPej5mPAxERGJgaEE+EpgspmNN7M04EpgSWTKEhGR/gy6heKcC5nZ54Bn6J1GeL9zbl3EKhMRkQ80pB64c+6PwB8jVIuIiBwF3UovIpKgFOAiIglKAS4ikqAU4CIiCcqcG9S9NYM7mVkDsH2Q314INEawnESg9zw86D0nv6G+32Occ0WHHoxpgA+Fma1yzlV5XUcs6T0PD3rPyS9a71ctFBGRBKUAFxFJUIkU4Iu8LsADes/Dg95z8ovK+02YHriIiPy9RBqBi4jIQRTgIiIJKiEC3MzOM7ONZrbJzG7zup5oMrOxZva8ma03s3Vm9kWva4oVM/Ob2Rtm9qTXtcSCmY0ws0fM7B0z22Bmp3hdU7SZ2ZfC/67fNrPfmFnSbWdlZvebWb2ZvX3QsQIz+7OZVYc/jozEueI+wA/aPPl8oBK4yswqva0qqkLAV5xzlcDJwE1J/n4P9kVgg9dFxNCPgaedc1OBmST5ezezcuALQJVzbjq9y1Bf6W1VUfFL4LxDjt0GLHXOTQaWhp8PWdwHOAdtnuyc6wT6Nk9OSs65Gufc6+HHzfT+T13ubVXRZ2ZjgI8C93pdSyyYWT7wYeA+AOdcp3Nuv6dFxUYKkGlmKUAWsNvjeiLOOfcCcOgu7AuAxeHHi4FLInGuRAjww22enPSBBmBmFcBsYIXHpcTCj4BbgB6P64iV8UAD8Itw2+heM8v2uqhocs7tAr4PvAfUAE3OuWe9rSpmSpxzNeHHtUBJJF40EQJ8WDKzHOD3wM3OuYDX9USTmV0I1DvnVntdSwylAHOAnznnZgOtROjX6ngV7vsuoPeHVxmQbWbXeFtV7LneudsRmb+dCAE+7DZPNrNUesP7Qefco17XEwOnAReb2TZ6W2TzzOwBb0uKup3ATudc329Xj9Ab6MnsHGCrc67BOdcFPAqc6nFNsVJnZqUA4Y/1kXjRRAjwYbV5spkZvX3RDc65H3pdTyw4577mnBvjnKug97/vX5xzST0yc87VAjvM7NjwobOB9R6WFAvvASebWVb43/nZJPmF24MsARaGHy8EnojEiw5pT8xYGIabJ58GXAu8ZWZrwsf+Jbz/qCSXzwMPhgcmW4BPe1xPVDnnVpjZI8Dr9M62eoMkvKXezH4DnAkUmtlO4HbgTuB3ZnYdvUtqXx6Rc+lWehGRxJQILRQRETkMBbiISIJSgIuIJCgFuIhIglKAi4gkKAW4iEiCUoCLiCSo/wPd5EH5JkKGgwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "array_of_radii = np.linspace(0,L,no_of_shells)\n",
    "dist_array = dist_array/(counter*len(dist_array))\n",
    "plt.plot(array_of_radii, dist_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "SiDbAwSl3n0o",
    "outputId": "52387c18-10d5-4f4a-d929-651594954d36",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "len(pair_distance_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4YwjKscR3n0p",
    "outputId": "e12a26ed-aea5-4aab-a670-b47190ca4412"
   },
   "outputs": [],
   "source": [
    "def create_picture(positions):\n",
    "    plt.cla()\n",
    "    plt.gca().set_aspect('equal')\n",
    "    plt.axis([0, L, 0, L])\n",
    "    plt.setp(plt.gca(), xticks=[0, L], yticks=[0, L])\n",
    "    for x,y in positions:\n",
    "        atom = plt.Circle((x, y), 0.45, fc='r')\n",
    "        plt.gca().add_patch(atom)\n",
    "create_picture(pos)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "9yJhSlWM3n0q"
   },
   "source": [
    "# Radial Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "4dcjvb893n0r"
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "4GlUXMT63n0s",
    "outputId": "f8b5c7f6-6f33-403a-cffc-067744c8e87b",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "radial_dis_data = []\n",
    "pair_distance_data = []\n",
    "for i in range(len(pos)):\n",
    "    for j in range(len(pos)):\n",
    "        if i!= j:\n",
    "            rij = pos[j]-pos[i]\n",
    "            pair_distance = np.dot(rij,rij)\n",
    "            pair_distance_data.append(pair_distance)\n",
    "            radial_dis_data.append(rdf(pair_distance))\n",
    "\n",
    "print(radial_dis_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 282
    },
    "id": "v5DcbodR6aKQ",
    "outputId": "0f7bbe6f-9955-4576-b949-84b2b609db7f"
   },
   "outputs": [],
   "source": [
    "plt.plot(pair_distance_data,radial_dis_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "qKLFlklV3n0s",
    "outputId": "f83f13f5-8db8-46b2-c9b0-2864865443c7"
   },
   "outputs": [],
   "source": [
    "pos, v = evolution_with_thermostate(pos, v)\n",
    "pos, v = evolution_without_thermostate(pos, v)\n",
    "data_for_rdf = []\n",
    "pos, v = evolution_without_thermostate(pos, v)\n",
    "data_for_rdf.append(rdf(pos))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "PE3Iys1M3n0s"
   },
   "outputs": [],
   "source": [
    "radiii = np.linspace(0, L, len(data_for_rdf))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "RQPFcsqC3n0t",
    "outputId": "2b490f0a-28bf-426b-db70-5120b2aba162"
   },
   "outputs": [],
   "source": [
    "data_for_rdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NBQ3NJRp3n0t"
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "radial distribution.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}