{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from ase.io import read\n",
    "from ase.visualize import view\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma   = 3.405\n",
    "epsilon = 119.8*8.616733e-5 # eV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def g_step(box, n_bin, d):\n",
    "    box    += 1\n",
    "    box    *= np.sqrt(2)/sigma\n",
    "    del_bin = box/n_bin  #width of a bin\n",
    "    \n",
    "    g    = np.zeros([n_bin,n_bin])\n",
    "    g[0] = np.linspace(0,box,n_bin)\n",
    "    g[0]+= del_bin*0.5\n",
    "    \n",
    "    for i in range(len(d)):\n",
    "        for j in range(i+1,len(d)):\n",
    "            g_index = d[i][j]/(del_bin)\n",
    "            g[1,g_index.astype(int)] += 2\n",
    "            # we count both for i and j.\n",
    "            \n",
    "    \n",
    "    #nomalized by number of particles\n",
    "    g[1] /= len(d[0])\n",
    "    #and by bin volume\n",
    "    \n",
    "    for k in range(len(g[1])):\n",
    "        g[1][k] /= ((k+1)**3 -k**3)*del_bin**3\n",
    "        \n",
    "    #side of the optimized cell: 5,269\n",
    "    vol = (5.269/sigma)**3\n",
    "    rho = 4./vol\n",
    "    g[1] /= np.pi*(4/3)*rho\n",
    "    return g;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 38 atom"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trajectory = read('production_38-pos-1.xyz', index='::2')\n",
    "N = len(trajectory[0])    #number of atoms\n",
    "n_step = len(trajectory)  #number of step\n",
    "\n",
    "#create a distance array\n",
    "dist = np.empty([0,N,N])\n",
    "for frame in trajectory:\n",
    "    dist = np.append(dist, [frame.get_all_distances()],axis=0)\n",
    "\n",
    "#Lennard Jones units\n",
    "\n",
    "dist  /= sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#understand what the distance array is:\n",
    "\n",
    "#dist0 = trajectory[0].get_all_distances()\n",
    "# print(dist0.shape)\n",
    "# print(dist0)\n",
    "# print(dist0[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_bin = 200  #number of bin in the graph\n",
    "box   = 25   #side of the bon in the simulation \n",
    "g_38  = np.zeros([n_bin,n_bin])\n",
    "\n",
    "#average during the thermalized part of the simulation\n",
    "for i in range(n_step):\n",
    "    g_38 += g_step(box,n_bin,dist[i])\n",
    "g_38 /= n_step\n",
    "\n",
    "plt.plot(g_38[0],g_38[1])\n",
    "plt.xlim(0,box/(2*sigma))\n",
    "plt.xlabel('distance [ units of $\\sigma$]')\n",
    "plt.ylabel('radial distribution function')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 150 atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trajectory = read('production_150-pos-1.xyz', index='::2')\n",
    "N = len(trajectory[0])    #number of atoms\n",
    "n_step = len(trajectory)  #number of step\n",
    "\n",
    "dist = np.empty([0,N,N])\n",
    "for frame in trajectory:\n",
    "    dist = np.append(dist, [frame.get_all_distances()],axis=0)\n",
    "dist  /= sigma\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_bin = 400  \n",
    "box   = 30   \n",
    "g_150 = np.zeros([n_bin,n_bin])\n",
    "\n",
    "\n",
    "for i in range(n_step):\n",
    "    g_150 += g_step(box,n_bin,dist[i])\n",
    "g_150 /= n_step\n",
    "\n",
    "plt.plot(g_150[0],g_150[1])\n",
    "plt.xlim(0,box/(2*sigma))\n",
    "plt.xlabel('distance [ units of $\\sigma$]')\n",
    "plt.ylabel('radial distribution function')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 450 atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trajectory = read('production_450-pos-1.xyz', index='::2')\n",
    "N = len(trajectory[0])    #number of atoms\n",
    "n_step = len(trajectory)  #number of step\n",
    "\n",
    "dist = np.empty([0,N,N])\n",
    "for frame in trajectory:\n",
    "    dist = np.append(dist, [frame.get_all_distances()],axis=0)\n",
    "dist  /= sigma\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_bin = 400  \n",
    "box   = 40   \n",
    "g_450 = np.zeros([n_bin,n_bin])\n",
    "\n",
    "\n",
    "for i in range(n_step):\n",
    "    g_450 += g_step(box,n_bin,dist[i])\n",
    "g_450 /= n_step\n",
    "\n",
    "plt.plot(g_450[0],g_450[1])\n",
    "plt.xlim(0,box/(2*sigma))\n",
    "plt.xlabel('distance [ units of $\\sigma$]')\n",
    "plt.ylabel('radial distribution function')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 817 atoms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trajectory = read('production_817-pos-1.xyz', index='::2')\n",
    "N = len(trajectory[0])    #number of atoms\n",
    "n_step = len(trajectory)  #number of step\n",
    "\n",
    "dist = np.empty([0,N,N])\n",
    "for frame in trajectory:\n",
    "    dist = np.append(dist, [frame.get_all_distances()],axis=0)\n",
    "dist  /= sigma\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_bin = 400  \n",
    "box   = 45   \n",
    "g_817 = np.zeros([n_bin,n_bin])\n",
    "\n",
    "\n",
    "for i in range(n_step):\n",
    "    g_817 += g_step(box,n_bin,dist[i])\n",
    "g_817 /= n_step\n",
    "\n",
    "plt.plot(g_817[0],g_817[1])\n",
    "plt.xlim(0,box/(2*sigma))\n",
    "plt.xlabel('distance [ units of $\\sigma$]')\n",
    "plt.ylabel('radial distribution function')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Thermalization plots (NVT)\n",
    "\n",
    "in the .energ file you can find\n",
    "\n",
    "- column 0: step\n",
    "- column 1: time \n",
    "- column 2: kinetic energy\n",
    "- column 3: temperature\n",
    "- column 4: potential energy\n",
    "\n",
    "Try to plot these values and describe what you see and what is conserved"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data38  = np.loadtxt('thermalization_38-1.ener' )\n",
    "data150 = np.loadtxt('thermalization_150-1.ener')\n",
    "data450 = np.loadtxt('thermalization_450-1.ener')\n",
    "data817 = np.loadtxt('thermalization_817-1.ener')\n",
    "\n",
    "data38  = np.transpose(data38)\n",
    "data150 = np.transpose(data150)\n",
    "data450 = np.transpose(data450)\n",
    "data817 = np.transpose(data817)\n",
    "\n",
    "plt.plot(data38[1],data38[3], label='38')\n",
    "plt.plot(data150[1],data150[3], label='150')\n",
    "plt.plot(data450[1],data450[3], label='450')\n",
    "plt.plot(data817[1],data817[3], label='817')\n",
    "\n",
    "plt.legend( loc='upper right')\n",
    "\n",
    "plt.xlabel(' time [fs]')\n",
    "plt.ylabel('T [K]')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Production plots (NVE)\n",
    "\n",
    "in the .energ file you can find\n",
    "\n",
    "- column 0: step\n",
    "- column 1: time \n",
    "- column 2: kinetic energy\n",
    "- column 3: temperature\n",
    "- column 4: potential energy\n",
    "\n",
    "Try to plot these values and describe what you see and what is conserved"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data38  = np.loadtxt('production_38-1.ener' )\n",
    "data150 = np.loadtxt('production_150-1.ener')\n",
    "data450 = np.loadtxt('production_450-1.ener')\n",
    "data817 = np.loadtxt('production_817-1.ener')\n",
    "\n",
    "data38  = np.transpose(data38)\n",
    "data150 = np.transpose(data150)\n",
    "data450 = np.transpose(data450)\n",
    "data817 = np.transpose(data817)\n",
    "\n",
    "plt.plot(data38[1],data38[3], label='38')\n",
    "plt.plot(data150[1],data150[3], label='150')\n",
    "plt.plot(data450[1],data450[3], label='450')\n",
    "plt.plot(data817[1],data817[3], label='817')\n",
    "\n",
    "plt.legend( loc='upper right')\n",
    "\n",
    "plt.xlabel(' time [fs]')\n",
    "plt.ylabel('T [K]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}