Part1_BarrierCrossing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Molecular Dynamics Asignment: Barrier Crossing In 1D\n",
    "Authors: *Stefan Hervø-Hansen*<sup>🐺</sup>, *Vidar Aspelin*<sup>🦡</sup>, and *Samuel Stenberg*<sup>🦄</sup>.  \n",
    "Division of Theoretical Chemistry, Lund University. \n",
    "\n",
    "Electronic Address:  \n",
    "🐺 stefan.hervo_hansen@teokem.lu.se  \n",
    "🦡 vidar.aspelin@teokem.lu.se  \n",
    "🦄 samuel.stenberg@teokem.lu.se\n",
    "\n",
    "<figure>\n",
    "    <img src=\"include/black_swan.png\" width=280 />\n",
    "</figure>\n",
    "<center>How does the title and asignment relate to the picture? The answer will be more aparent by the end of the asignment.</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Intended learning outcomes\n",
    "1. The effects on the time step on the total energy fluctuations\n",
    "2. Transformation of energy (potential vs kinetic energy)\n",
    "3. Concept of phase space and available microstates in different ensembles\n",
    "4. Sampling differences between MD and MC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theory and Questions for Reflection\n",
    "Consider the movement of a single particle that moves on a 1D potential energy surface with the following functional form:\n",
    "$$\n",
    "U(x)=\n",
    "\\begin{cases}\n",
    "\\epsilon B x^2, & \\text{if}\\ x \\leq 0.0\\\\\n",
    "\\epsilon x_{\\circ}^2 \\left[ 1 - \\cos \\left( \\frac{2\\pi x}{x_\\circ} \\right)\\right], & \\text{if}\\ 0.0 \\leq x \\leq 1.0\\\\\n",
    "\\epsilon B (x_{\\circ}-x)^2, & \\text{if}\\ x \\geq 1.0\n",
    "\\end{cases}\n",
    "\\tag{eq. 1}\n",
    "$$\n",
    "The energy, force and the derivative of the force should be continuous functions of the position $x$ with $\\epsilon$ > 0 and $x_\\circ$ is a standard position set to 1 for units to cancel, here $x_\\circ=1 \\ \\mathrm{nm}$. We want to investigate this system using Molecular dynamics.\n",
    "\n",
    "**Question 1: Briefly explain Molecular dynamics (include the concepts Newton's equations of motion, trajectory, integrator, time step, NVE ensemble, energy)**  \n",
    "\n",
    "In order to propagate positions and velocities in time, we need information about the acceleration, which is obtained from the forces. To make sure that our system is stable and can be propagated correctly, we need to make the potential, the force, and the derivative of the force (with respect to x) continuous.\n",
    "\n",
    "<figure>\n",
    "<img src=\"include/MC_vs_MD_potential.png\" width=800 />\n",
    "</figure>\n",
    "<center>A comparison between a non-continuous potential (hard spheres) and a continuous potential (Weeks-Chandler-Andersen.)</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__Question 2: Derive an expression for *B* so that the potential, force and derivative of the force becomes continuous and plot/sketch the energy landscape. Also insert your derived expression for *B* in the code containing the function Potential below.__  \n",
    "\n",
    "The program pymd is provided that integrates the equation of motion of the particle on the potential energy surface, starting at the position 0 ($x(t=0)=0$), using either molecular dynamics and methods thereof or Monte Carlo simulation. The program functionalities has been listed below:  \n",
    "<ul style=\"list-style-type:none;\">\n",
    "  <li>Molecular dynamics, NVE ensemble</li>\n",
    "  <li>Molecular dynamics with Andersen thermostat</li>\n",
    "  <li>Molecular dynamics with Berendsen thermostat</li>\n",
    "  <li>Monte Carlo </li>\n",
    "</ul>  \n",
    "For the functionalities with thermostats, we are approximating the NVT ensemble using two different thermostats, in particular the Andersen thermostat and Berendsen thermostat which can be categorized as a stochastic method and velocity rescale method respectively.\n",
    "\n",
    "\n",
    "Before conducting simulations, review the part of the code (pymd.py) containing the Velocity Verlet integrator by executing the cell below, and try to understand what is happening."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!sed -n 34,39p pymd.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import random\n",
    "import pymd as sim # This is where we import the program used in this laboratory session"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function defining the potential (eq.1) and force (eq.2) \n",
    "eps = 15\n",
    "B = # FINISH THE CODE: Add expression for B\n",
    "def Potential(x):\n",
    "    '''Function returning energy and force depending on the position, x, on the 1D energy surface.'''\n",
    "    if x <= 0.0:\n",
    "        U = eps *  B * x**2\n",
    "        F = eps * -2 * B * x\n",
    "    elif x >= 1.0:\n",
    "        U = eps * B * (x-1.0)**2\n",
    "        F = eps * -2 * B * (x-1.0)\n",
    "    else:\n",
    "        U = eps * (1.0 - np.cos(2*np.pi*x))\n",
    "        F = eps * -2.0 * np.pi * np.sin(2.0*np.pi*x)\n",
    "    return U, F      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the energy and force here.\n",
    "# Plot energy:\n",
    "x = np.arange(-1, 2, 0.01)\n",
    "y = [Potential(x_)[0] for x_ in x]\n",
    "plt.plot(x, y)\n",
    "plt.xlabel('$x$ (nm)')\n",
    "plt.ylabel('$U$(x) (kJ/mol)')\n",
    "\n",
    "# FINISH THE CODE: Plot the force. Check so that the force is everywhere differentiable."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question 3: Conduct short simulations of the NVE ensemble using different time steps but the same simulation time and comment on the differences in energy profiles (potential, kinetic and total energy). Elaborate on the differences you observe!**\n",
    "\n",
    "Newton's equations tell us that if we have a set of particle positions and particle velocities at an arbitrary time, it is possible to predict the positions and velocities at any one particular instant in time. Therefore, classical dynamics of a $N$-particle system can be expressed by specifying the full set of 6$N$ variables ($2N$ for each dimension). The $6N$ variables can be regarded as a single point in a $6N$-dimensional space named the _phase space_. The phase space is a Cartesian space and solving Newton's equations yields a trajectory in phase space. Consequently, classical motion can be described by the motion of a point along a trajectory in phase space, with each point in phase space representing a possible microstate. \n",
    "\n",
    "To be able to visualize phase space, we will in the following plot the phase space for 1 particle in 1 dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SIMULATION SETTINGS, MD-simulations (NVE)\n",
    "Temperature = 300                                 # [K]\n",
    "Nparticles = 100                                  # Number of independent particles\n",
    "Tstep = 0.001                                      # [ps]\n",
    "SimTime = 1                                       # Simulation time [ps]\n",
    "Nstep = int(SimTime/Tstep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nve = sim.Simulator(Nstep, Nparticles, Tstep, Temperature)\n",
    "\n",
    "for j in range(nve.Nparticles): # Particle loop  \n",
    "    nve.Velocity[j] = 1 # md.RandomVelocity(Temperature, m)\n",
    "\n",
    "\n",
    "for i in range(Nstep): # Main Loop\n",
    "    nve.VelocityVerlet_NVE(Potential) # integrate motion of particles\n",
    "    nve.SampleData() # sample data for step i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting potential, kinetic and total energy as a function of simulation time\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "\n",
    "ax.plot(nve.Data_time[:], np.sum(nve.Data_Epot[:], axis=1), label=r'Potential energy', lw=3, alpha=0.7)\n",
    "ax.plot(nve.Data_time[:], np.sum(nve.Data_Ekin[:], axis=1), label=r'Kinetic energy', lw=3, alpha=0.7)\n",
    "ax.plot(nve.Data_time[:], np.sum(nve.Data_Etot[:], axis=1), label=r'Total energy', lw=3, alpha=0.7)\n",
    "\n",
    "# Graphics settings\n",
    "ax.minorticks_on()\n",
    "ax.tick_params(axis='both',which='minor',length=5,width=2,labelsize=18)\n",
    "ax.tick_params(axis='both',which='major',length=8,width=2,labelsize=18)\n",
    "ax.legend(fontsize=16, bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "\n",
    "ax.set_title(r'Molecular Dynamics Stability', fontsize=22)\n",
    "ax.set_xlabel(r'Time (ps)', fontsize=22, labelpad=20)\n",
    "ax.set_ylabel(r'Energy (kJ/mol)', fontsize=22, labelpad=20)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question 4: What do you expect the phase space trajectories for the 1 particle to look like in terms of density and shape for the following situations, given a (i) Harmonic potential and (ii) the potential given in Eq. 1? _Hint_: The first term in Eq. 1 is a harmonic oscillator. How would the shape of phase space change if we added more energy to the NVE ensemble?**\n",
    "1. *NVE* ensemble.\n",
    "2. *NVT* ensemble."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question 5: Using the simulations methods (NVE, NVT and MC) for low temperatures (150 < T < 300 K), comment on the following:**\n",
    "\n",
    "*Note*: When simulating in MC, we no longer have velocities. This means that one dimension in phase space, the velocity dimension, dissapears and we are left with only the positions (the *configurational* dimension). So when talking about phase space in MC, we say configurational space rather than phase space. For the *NVT* ensemble, we will only be using the Andersen thermostat in this lab exercise.\n",
    "\n",
    "1. Does the result correspond with what you expected in terms of the density and shape? Comment on any differences. \n",
    "2. Why does the phase space distribution for the *NVE* and *NVT* ensemble look so different?\n",
    "3. Why does the configurational space distribution of the MC scheme look so much different compared to the phase space distribution of the MD schemes at reasonably low temperatures (<10000K)?\n",
    "4. See cell with question 5.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A look into pymd for the Andersen thermostat random velocity asignment\n",
    "!sed -n 84,89p pymd.py\n",
    "!echo \" \"\n",
    "!sed -n 110,112p pymd.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SIMULATION SETTINGS, MD-simulations (NVE and NVT)\n",
    "Temperature = 300                                 # [K]\n",
    "Nparticles = 10\n",
    "# Number of independent particles\n",
    "Tstep = 0.001                                      # [ps]\n",
    "SimTime = 500                                       # Simulation time [ps]\n",
    "Nstep = int(SimTime/Tstep)                        # Calculated number of simulation steps \n",
    "MU = 1                                            # Collision frequency (used for NVT-ensemble simulation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MD-simulations with NVE and NVT\n",
    "\n",
    "nve = sim.Simulator(Nstep, Nparticles, Tstep, Temperature)\n",
    "nvt = sim.Simulator(Nstep, Nparticles, Tstep, Temperature, MU)\n",
    "\n",
    "for j in range(nve.Nparticles): # Particle loop  \n",
    "    nve.Velocity[j] = nve.RandomVelocity(Temperature, 1) #Set random velocity from MB distribution\n",
    "    nvt.Velocity[j] = nvt.RandomVelocity(Temperature, 1)\n",
    "\n",
    "for i in range(Nstep): # Main Loop for NVE simulation\n",
    "    nve.VelocityVerlet_NVE(Potential) # integrate motion of particles\n",
    "    nve.SampleData() # sample data for step i\n",
    "    \n",
    "for i in range(Nstep): # Main Loop for NVT simulation\n",
    "    nvt.VelocityVerlet_NVE(Potential) # integrate motion of particle j\n",
    "    nvt.Thermostat_And() # Applying thermostat to simulate in NVT\n",
    "    nvt.SampleData() # sample data for step i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting phase space distributions for MD simulations (NVE and NVT)\n",
    "fig, (ax1, ax2) = plt.subplots(ncols = 2, figsize=(12,6))\n",
    "\n",
    "ax1.hist2d(nve.Data_pos.ravel(), nve.Data_vel.ravel(), density=True, bins=100)    \n",
    "    \n",
    "# Graphics settings\n",
    "ax1.minorticks_on()\n",
    "ax1.tick_params(axis='both',which='minor',length=5,width=2,labelsize=18)\n",
    "ax1.tick_params(axis='both',which='major',length=8,width=2,labelsize=18)\n",
    "\n",
    "ax1.set_title(r'Phase Space Distribution, NVE', fontsize=22)\n",
    "ax1.set_ylabel(r'Velocity (nm/ps)', fontsize=22, labelpad=20)\n",
    "ax1.set_xlabel(r'Position (nm)', fontsize=22, labelpad=20)\n",
    "\n",
    "ax2.hist2d(nvt.Data_pos.ravel(), nvt.Data_vel.ravel(), density=True, bins=100)    \n",
    "    \n",
    "# Graphics settings\n",
    "ax2.minorticks_on()\n",
    "ax2.tick_params(axis='both',which='minor',length=5,width=2,labelsize=18)\n",
    "ax2.tick_params(axis='both',which='major',length=8,width=2,labelsize=18)\n",
    "\n",
    "ax2.set_title(r'Phase Space Distribution, NVT', fontsize=22)\n",
    "ax2.set_ylabel(r'Velocity (nm/ps)', fontsize=22, labelpad=20)\n",
    "ax2.set_xlabel(r'Position (nm)', fontsize=22, labelpad=20)\n",
    "\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SIMULATION SETTINGS, MC-simulation\n",
    "Temperature = 300                                 # [K]\n",
    "Nparticles = 10                                  # Number of independent particles\n",
    "Tstep = 0.001                                     # [ps]\n",
    "SimTime = 500                                       # Simulation time [ps]\n",
    "Nstep = int(SimTime/Tstep)                        # Calculated number of simulation steps \n",
    "dp = 1                                            # Displacement parameter for MC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MC-simulations\n",
    "\n",
    "mc = sim.Simulator(Nstep, Nparticles, Tstep, Temperature, dp=dp)\n",
    "\n",
    "for i in range(Nstep): # Main Loop for NVE simulation\n",
    "    for j in range(mc.Nparticles): # Particle loop  \n",
    "        mc.MonteCarlo(j, Potential) # integrate motion of particle j\n",
    "    mc.SampleData() # sample data for step i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting phase configurational space distribution for MC-simulation\n",
    "fig, ax = plt.subplots(figsize=(7,6))\n",
    "\n",
    "ax.hist(mc.Data_pos.ravel(), density=True, bins=100)\n",
    "\n",
    "\n",
    "    \n",
    "# Graphics settings\n",
    "ax.minorticks_on()\n",
    "ax.tick_params(axis='both',which='minor',length=5,width=2,labelsize=18)\n",
    "ax.tick_params(axis='both',which='major',length=8,width=2,labelsize=18)\n",
    "\n",
    "ax.set_title(r'Configurational Space Distribution',fontsize=22)\n",
    "ax.set_ylabel(r'Probability', fontsize=22, labelpad=20)\n",
    "ax.set_xlabel(r'Position (nm)', fontsize=22, labelpad=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**5.4.** Show that the MD data obtained using the Andersen thermostat matches the analytical probability distribution for velocities and positions for the harmonic oscillator. Use the code provided below, and insert the missing expression for the partition function. *Hint*:\n",
    "$$\n",
    "\\int_{-\\infty}^{\\infty} e^{-a(x+b)^2} \\mathrm{d}x = \\sqrt{\\frac{\\pi}{a}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "################################################\n",
    "#                QUESTION 5.4                  #\n",
    "# Boltzman distribution: Potential energy part #\n",
    "################################################\n",
    "\n",
    "plt.hist(nvt.Data_pos.ravel(), density=True, bins=1000)\n",
    "\n",
    "# FINISH THE CODE: START\n",
    "# Potential energy partition function:\n",
    "partition_function_PE =  #insert partition function here\n",
    "\n",
    "def probability_PE(m, Temperature, x):  # Probability as a function of position \n",
    "    return #insert probability here\n",
    "\n",
    "# FINISH THE CODE: END\n",
    "\n",
    "x = np.arange(-0.2, 0.2, 0.0001)\n",
    "plt.plot(x, probability_PE(1, Temperature, x), 'r--', linewidth=3)\n",
    "plt.xlabel(r'$\\bf{x}\\ \\mathrm{(nm)}$')\n",
    "plt.ylabel(r'$\\bf{P(x)}$')\n",
    "plt.show()\n",
    "# Self-Check: The sum of all probabilities should always equal 1!\n",
    "print('The area under the curve is equal to: {:03.2f}'.format(np.trapz(probability_PE(1,Temperature,x), x)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################################\n",
    "#                QUESTION 5.4                #\n",
    "# Boltzman distribution: Kinetic energy part #\n",
    "##############################################\n",
    "plt.hist(nvt.Data_vel.ravel(), density=True, bins=100)\n",
    "\n",
    "# FINISH THE CODE: START\n",
    "# Kinetic partition function:\n",
    "partition_function_KE =  #insert partition function here\n",
    "\n",
    "# Probability as a function of velocity\n",
    "def probability_KE(m, Temperature, v):\n",
    "    return  #insert probability here\n",
    "\n",
    "# FINISH THE CODE: END\n",
    "\n",
    "v = np.arange(-5, 5, 0.001)\n",
    "plt.plot(v, probability_KE(1, Temperature, v), 'r--', linewidth=3)\n",
    "plt.xlabel(r'$\\bf{v}\\ \\mathrm{(nm/ps)}$')\n",
    "plt.ylabel(r'$\\bf{P(v)}$')\n",
    "plt.show()\n",
    "\n",
    "# Self-Check: Should always equal 1!\n",
    "print('The area under the curve is equal to: {:3.2f}'.format(np.trapz(probability_KE(1,Temperature,v), v)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question 6: Calculate at which temperature half the particles would be able to cross the energy barrier and test your hypothesis using Molecular dynamics in the NVE ensemble. Comment on the result and compare with the mean speed of the Maxwell-Boltzmann distribution at the predicted temperature.**  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##############################################\n",
    "#                QUESTION 6                  #\n",
    "#         Prediction of temperature          #\n",
    "##############################################\n",
    "\n",
    "L = # FINISH THE CODE: Insert lower integrating limit\n",
    "\n",
    "T_list = np.arange(300,50000,100) # List of temperature\n",
    "v_bot = np.arange(0, 10000, 0.1)  # Limits for calculating total Boltzmann weigth\n",
    "v_top = np.arange(L, 10000, 0.1) # Limits for calculating the Boltzmann weigth of crossing\n",
    "frac_list = []\n",
    "\n",
    "# Integrand\n",
    "def integrand(Temperature, v):\n",
    "    return np.exp(-0.5*1*v**2/(nve.kB*Temperature))  # Kinetic boltzmann weigth\n",
    "\n",
    "# Integrate using trapezoidal rule\n",
    "for T in T_list:\n",
    "    p_cross = np.trapz(integrand(T,v_top), v_top) # Integrate to get Boltzmann weigth of crossing\n",
    "    p_tot = np.trapz(integrand(T,v_bot), v_bot)   # Integrate to get total boltzmann weigth\n",
    "    frac = p_cross / p_tot                        # Fraction\n",
    "    # Find temperature where fraction of the integrals is 0.5 i.e probability of half crossing\n",
    "    if (abs(frac-0.5)) < 0.0008:\n",
    "        print(r'The temperature at which half the particles can cross the barrier is {}'.format(T))\n",
    "    frac_list.append(frac)\n",
    "    \n",
    "plt.plot(T_list, frac_list)\n",
    "\n",
    "plt.xlabel(\"Temperature (Kelvin)\")\n",
    "plt.ylabel(\"Fraction of particles\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SIMULATION SETTINGS, MD-simulations (NVE)\n",
    "Temperature =                                     # FINISH THE CODE: Insert calculated temperature\n",
    "Nparticles = 10000                                # Number of independent particles\n",
    "Tstep = 0.001                                     # [ps]\n",
    "SimTime = 1                                       # Simulation time [ps]\n",
    "Nstep = int(SimTime/Tstep)                        # Number of simulation steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MD-simulations with NVT\n",
    "\n",
    "nve = sim.Simulator(Nstep, Nparticles, Tstep, Temperature)\n",
    "\n",
    "for j in range(nve.Nparticles):\n",
    "    nve.Velocity[j] = nve.RandomVelocity(Temperature, 1)  # Set random velocity\n",
    "    \n",
    "for i in range(Nstep): # Main Loop for NVT simulation\n",
    "    nve.VelocityVerlet_NVE(Potential) # integrate motion of particles\n",
    "    nve.SampleData() # sample data for step i\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "particle_id = np.arange(0, Nparticles, 1)\n",
    "particle_id = np.stack((particle_id, [False]*Nparticles), axis=-1)\n",
    "Tcounter = 0\n",
    "\n",
    "# Check which particles has crossed the barrier\n",
    "for particle in range(Nparticles):\n",
    "    for position in nve.Data_pos.T[particle]:\n",
    "        if position > 0.5:\n",
    "            particle_id[particle][1] = True\n",
    "            break\n",
    "\n",
    "# Count the number of particles which has crossed the barrier\n",
    "for particle in particle_id:\n",
    "    if particle[1] == True:\n",
    "        Tcounter += 1\n",
    "\n",
    "print('{}% of the particles crossed the barrier at T = {} K'.format(Tcounter/Nparticles*100, Temperature))"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}