Part2_Diffusion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diffusion from displacement of particles\n",
    "In this part of the lab we are going to calculate the self-diffusion of water from data of the positions over time. In the program, a textfile called _msd.txt_ is generated after a simulation has been run, containing the mean squared displacement of water as a function of time. The mean squared displacement of water is calculated by summing the squares of the oxygen displacements from their respective reference positions determined by their start positions, $x_{0,n}$, at predefined points in time. Then, for each sampled point in time, the sum is divided by the total number of water molecules,\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle(x(t)-x(t=0))^2\\rangle = \\frac{1}{N}\\sum\\limits_{n=1}^N (x_n(t)-x_n(t=0))^2\n",
    "\\label{msd}\n",
    "\\tag{1}\n",
    "\\end{align}\n",
    "$$\n",
    "where $N$ is the total number of particles. \n",
    "\n",
    "If particles interact, they will establish a random walk behaovior after enough \"collisions\" have ocuured, causing the mean squared displacement to increase lineary with time. However, in the start of a simulation, it is common that one observes a non-linear dependence of the mean squared displacement on time. This is known as the _ballistic_ regime, taking place because not enough collisions have occured to establish the random walk behavior. Mathematically, this can be shown by considering two different cases,\n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "x_1(t+\\Delta t) = x(t) \\pm \\delta\n",
    "\\label{diff1}\n",
    "\\tag{2}\n",
    "\\end{equation}\n",
    "$$   \n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "x_2(t+\\Delta t) = x(t) + \\delta\n",
    "\\label{diff2}\n",
    "\\tag{3}\n",
    "\\end{equation}\n",
    "$$   \n",
    "where equation \\ref{diff1} describes random walk motion, with equal probabilities to move particles in the positive and negative direction, whereas equation \\ref{diff2} describes ballistic motion, only displacing particles in the positive direction. The displacement parameter, $\\delta$, is assumed to keep constant over time.\n",
    "\n",
    "2. Run simulations with 216 water molecules for a series of temperatures (4-5 different, you choose which) and determine the diffusion coefficient of water for each temperature using the generated textfile _msd.txt_. Plot the mean squared displacements for all temperatures in the same plot, calculate the corresponding diffusion coefficients and report the values in a table. \n",
    "\n",
    "*Hints*:\n",
    "\n",
    "i) The diffusion coeficient can be obtained using the relation,\n",
    "$$\n",
    "\\begin{equation}\n",
    "D = \\frac{\\langle(x(t)-x_0)^2\\rangle}{2 Nt}\n",
    "\\tag{4}\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "ii) The curve for mean squared displacement as a function of time is often not linear. To obtain a reasonable value of the diffusion coefficient, one most commonly makes a linear regression of the part of the curve that shows nearly linear behavior. This is done below using the package stats of scipy:\n",
    "\n",
    "```python\n",
    "slope, intercept, r_value, p_value, std_err = stats.linregress(time_array,msd_array)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading packages\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import shutil\n",
    "import os\n",
    "from scipy import stats\n",
    "from tempfile import mkstemp\n",
    "from os import fdopen, remove\n",
    "# Importing modules for loading and viewing\n",
    "import mdtraj as md\n",
    "import nglview as ngl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Naming directories to be used\n",
    "datadir = 'data/'\n",
    "!pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining function to replace the value of input parameters in the input file 'config.mo' \n",
    "def replace(file_path, pattern, subst):\n",
    "    with open(file_path,\"r+\") as f:\n",
    "        new_f = f.readlines()\n",
    "        f.seek(0)\n",
    "        for line in new_f:\n",
    "            if pattern not in line:\n",
    "                f.write(line)\n",
    "            else:\n",
    "                f.write(subst+\"\\n\")\n",
    "        f.truncate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Ts = [ ] # ENTER CODE HERE: Enter the temperatures you want to simulate, separated by a comma\n",
    "iterations =  # multiply this with the timestep (0.001 ps) and you get the simulation time in ps. Use iterations 5000 <= x <= 50000 \n",
    "\n",
    "replace('data/config.mo', 'iter', 'iter '+str(iterations))\n",
    "replace('data/config.mo', 'simple', 'simple '+str(0))\n",
    "replace('data/config.mo', 'timestep',    'timestep '    + str(0.001))\n",
    "\n",
    "for T in Ts:\n",
    "    print(T)\n",
    "    replace('data/config.mo', 'temperature', 'temperature '+str(T))\n",
    "    !./md/src/main data/WAT_216.equil.gro data/WAT_216.par data/config.mo\n",
    "    os.rename('msd.txt', 'msd_T_'+str(T)+'.txt')\n",
    "    shutil.move(\"msd_T_\"+str(T)+\".txt\", datadir+\"msd_T_\"+str(T)+\".txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for T in Ts: # Looping over temperatures\n",
    "    time = np.loadtxt(datadir+\"msd_T_\"+str(T)+\".txt\", usecols=0) # Loading time array for temperature T\n",
    "    msd = np.loadtxt(datadir+\"msd_T_\"+str(T)+\".txt\", usecols=1) # Loading mean squared displacement for each time, averaged over all oxygens, for temperature T\n",
    "\n",
    "    ### Plot the msd vs time\n",
    "    plt.plot() # ENTER CODE HERE: enter arguments to plot and add labels with label=\"arbitraryLabelName\"\n",
    "plt.xlabel('$time$ / ps', fontsize=13) # Assigning x-axis label\n",
    "plt.ylabel(r'$\\langle(x(t)-x_0)^2\\rangle$ / (nm$^2$ ps$^{-1})$', fontsize=13) # Assigning y-axis label\n",
    "plt.legend() # Plotting legend\n",
    "plt.savefig(datadir+\"msd_vs_time.pdf\") # Saving figure as pdf in subdirectory \"data\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Plot the diffusion coefficient as a function of temperature in units cm$^2$/s and K, respectively. Discuss what you observe. Why does diffusion increase with temperature? Use keywords like entropy/thermal energy, enthalpy and local minima.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting diffusion vs temperature\n",
    "i = 0 # Initializing counter\n",
    "D = np.zeros([len(Ts)]) # Initializing Diffusion vector\n",
    "N =    # ENTER CODE HERE: Number of dimensions (to be used in the Einstein relation)\n",
    "\n",
    "for T in Ts: # Looping over temperatures\n",
    "    time = np.loadtxt(datadir+\"msd_T_\"+str(T)+\".txt\", usecols=0) # Loading time array for temperature T\n",
    "    msd = np.loadtxt(datadir+\"msd_T_\"+str(T)+\".txt\", usecols=1) # Loading mean squared displacement for each time, averaged over all oxygens, for temperature T\n",
    "    slope, intercept, r_value, p_value, std_err = stats.linregress(time,msd)\n",
    "    #print(msd)\n",
    "    D[i] =     # ENTER CODE HERE: enter the expression for diffusion according to eq. 4, using the slope calculated. You will need a conversion factor to convert the slope in nm^2/ps to cm^2/s. \n",
    "    print(\"D(T=\"+str(T)+\" K) = \", D[i])\n",
    "    i = i+1 # Increasing counter by one\n",
    "plt.plot()      # ENTER CODE HERE: Enter arguments (temperature and diffusion)\n",
    "plt.xlabel('$T$ / K')\n",
    "plt.ylabel('$D$ / (cm$^2$ s$^{-1}$)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualize\n",
    "# Modifying gro files to be compatible with nglview (takes some time: a couple of minutes...)\n",
    "for t in np.arange(0,10,1):\n",
    "    replace('output.gro', 'Generated by Slaymulator, t='+str(t), 'Generated by Slaymulator, t= '+str(t)+'.0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading trajectory\n",
    "waterTraj = md.load(\"output.gro\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading trajectory for viewing\n",
    "viewTraj = ngl.show_mdtraj(waterTraj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewTraj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Errors in MD Simulations:\n",
    "A computer can normally only handle 15 significant numbers. This results in a small error each time a \n",
    "calculation is performed. To get more insight into this phenomenon and it's implications we take advantage \n",
    "of the deterministic nature of MD. \n",
    "\n",
    "In this part you will run a short simulation, approx 5ps, of 12 ions initially placed in the xy-plane, using a structure file called \"2d.gro\". You should then use the last frame of this simulation to run another simulation but this time backwards, with the same number of steps as the first. Since MD is deterministic you should end up at exactly the same configuration as the one you started with. It is important that you choose the same number of steps in both simulations to be able to get a clear picture of the error and the error propagation.\n",
    "\n",
    "Since the ions will only move in two dimensions you are able to plot the positions of one particle in a graph. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run forward simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iterations = 170000 # multiply this with the timestep (0.001 ps) and you get the simulation time in ps\n",
    "\n",
    "# Modifying config.mo to contain what we want\n",
    "replace('data/config.mo', 'iter',        'iter '        + str(iterations))\n",
    "replace('data/config.mo', 'temperature', 'temperature ' + str(300))\n",
    "replace('data/config.mo', 'timestep',    'timestep '    + str(0.001))\n",
    "replace('data/config.mo', 'simple',        'simple '    + str(1))\n",
    "\n",
    "# Running\n",
    "!./md/src/main data/2D_NaCl.gro data/2D_NaCl.par data/config.mo\n",
    "\n",
    "# Renaming output files\n",
    "!mv output.gro forward.gro \n",
    "!mv track.txt track_forw.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run backward simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "replace('data/config.mo', 'timestep',    'timestep '    + str(-0.001)) # Editing timestep in config.mo to run backwards \n",
    "\n",
    "!tail -n 15 forward.gro | tee forwardLastFrame.gro # Copying last frame from forward.gro to forwardLastFrame.gro\n",
    "\n",
    "# Running\n",
    "!./md/src/main forwardLastFrame.gro data/2D_NaCl.par data/config.mo\n",
    "\n",
    "# Renaming output files\n",
    "!mv output.gro backward.gro \n",
    "!mv track.txt track_back.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading track files\n",
    "forward = np.loadtxt('track_forw.txt', unpack=True)\n",
    "backward = np.loadtxt('track_back.txt', unpack=True)\n",
    "\n",
    "fx = forward[1][1:] - forward[1][:-1]\n",
    "fy = forward[2][1:] - forward[2][:-1]\n",
    "gx = backward[1][1:] - backward[1][:-1]\n",
    "gy = backward[2][1:] - backward[2][:-1]\n",
    "\n",
    "# Plotting forward and backward trajectory\n",
    "plt.figure()\n",
    "plt.plot(forward[1], forward[2], label='forward')\n",
    "plt.plot(backward[1], backward[2], label='backward')\n",
    "plt.xlabel(\"x\")\n",
    "plt.xlabel(\"y\")\n",
    "plt.legend()\n",
    "\n",
    "# Plotting forward and backward trajectory\n",
    "#print(fx)\n",
    "print(forward[1][:-1])\n",
    "plt.figure()\n",
    "plt.plot(fy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Plot the positions of one particle from both simulations in the same graph. \n",
    "\n",
    "2. Compare the structure of the first frame in the first simulation and the last frame in the second simulation. \n",
    "\n",
    "3. Shortly discuss the results you have found."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(forward[1], forward[2], label='forward')\n",
    "plt.plot(backward[1], backward[2], label='backward')\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Modifying gro files to be compatible with nglview\n",
    "for t in np.arange(0,1000,1):\n",
    "    replace('backward.gro', 'Generated by Slaymulator, t='+str(t), 'Generated by Slaymulator, t= '+str(t)+'.0')\n",
    "    replace('forward.gro', 'Generated by Slaymulator, t='+str(t), 'Generated by Slaymulator, t= '+str(t)+'.0')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "forwardTraj = md.load(\"forward.gro\")\n",
    "backwardTraj = md.load(\"backward.gro\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading files for viewing\n",
    "viewForward = ngl.show_mdtraj(forwardTraj)\n",
    "viewBackward = ngl.show_mdtraj(backwardTraj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Viewing forward trajectory\n",
    "viewForward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Viewing backward trajectory\n",
    "viewBackward"
   ]
  },
  {
   "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}