{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ef0fcfdd",
   "metadata": {},
   "source": [
    "# Week XII - Trajectories: Dynamics Visualization and Radial Distribution Function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ba562f9",
   "metadata": {},
   "source": [
    "In this tutorial we analyze **trajectories** obtained by Monte Carlo (MC) simulations, in particular we focus on the visualization of the dynamics and the calculation of the radial distribution function (RDF). We will consider MC simulations of 2D hard disks, but the same considerations applies to molecular dynamics (MD) and other systems.\n",
    "\n",
    "> We will analyse the trajectories with a post-processing approach:\n",
    "> * The Fortran code `hd-MC.f90` will be use to run the dynamics. \n",
    "> * We have slightly modified `hd-MC.f90` to store the trajectory in a standardized format (ext XYZ)\n",
    "> * We will use Python to interface with existing visualization libraries and calculate the RDF on top of the stored `.xyz` trajectories"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23818a66",
   "metadata": {},
   "source": [
    "NB: we will use the Python package `ase` for parsing and visualization (you can install it with by ``pip install ase``)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7febb285",
   "metadata": {},
   "source": [
    "## Visualizaing the trajectory\n",
    "\n",
    "* MC and MD simulations typically involve the calculation of a number of quantities for analysis, such as the RDF, the mean square displacement (MSD), the autocorrelation function and more.\n",
    "* While these properties can be calculated \"on the fly\" during the simulations, it is often more convenient to store the trajectories and calculate the desired quantitites later. That mirrors exactly the situation in laboratory experiments: **data acquisition and data analysis are typically performed at different times**.\n",
    "* The same data analysis routines can be used to study different systems, potentially obtained with different simulations codes, as long as a standardized *data format* is used.\n",
    "* Hence, it is often very powerful to produce outputs in widely recognized formats and extensions, that allows to use existing packages for data analysis: you (or someone else) might have already produced the relevant function or subroutine."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bc6e34e",
   "metadata": {},
   "source": [
    "### Extended XYZ and ASE\n",
    "- We will use the extend XYZ format for storing the trajectories (see for instance https://www.ovito.org/docs/dev/reference/file_formats/input/xyz.html)\n",
    "- We will use the Python package Atomic Simulation Environment (ASE) for reading the XYZ file (https://wiki.fysik.dtu.dk/ase/) and visualize the dynamics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16e0f0d9",
   "metadata": {},
   "source": [
    "### Simulations and post processing\n",
    "1. First we look the code ``hd-MC_v2.f90``, which has been modified to store trajectories in the extended XYZ format (now we move to the code editor window)\n",
    "2. Let's run in the terminal some MC with the Fortran code `hd-MC_v2.x` (N=30, 6x5, 500 steps for equilibration, 5000 for sampling)\n",
    "3. Now we analyze the trajectoy ``traj.xyz``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79fc8618",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from ase.io import read\n",
    "# We load and parse the trajectory\n",
    "traj = read('./Fortran_code/traj.xyz',index=':')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96319a71",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ASE has converted the ext XYZ format in an internal format\n",
    "# Each element of the list is a snapshot of the trajectory\n",
    "traj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8708bf62",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = traj[0]\n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08ec488a",
   "metadata": {},
   "outputs": [],
   "source": [
    "_ = [a.set_pbc([True,True,False]) for a in traj]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7edc5d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ase.visualize import view\n",
    "# Here we visualize the trajectory\n",
    "view(traj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0e6592e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from ase.visualize.plot import plot_atoms\n",
    "import matplotlib.pyplot as plt\n",
    "fig, ax = plt.subplots()\n",
    "sn = traj[500]\n",
    "plot_atoms(sn,ax,radii=[0.5]*sn.get_global_number_of_atoms())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83777c5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#sn = traj[501]\n",
    "#plot_atoms(sn,ax,colors=['orange']*sn.get_global_number_of_atoms(),radii=[0.5]*sn.get_global_number_of_atoms())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb284632",
   "metadata": {},
   "source": [
    "## Radial Distribution Function (RDF)\n",
    "\n",
    "* The RDF in 2D counts the number of particle pairs which are separated by a distance between $r$ and $r+dr$\n",
    "\n",
    "$$g(r) = \\frac{2}{2\\pi r \\Delta r \\rho N}\\langle \\sum_{i=1}^{N-1}\\sum_{j>i}\\delta(r-r_{ij}) \\rangle$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48326716",
   "metadata": {},
   "source": [
    "> Now we implement the RDF calculation on top of the stored trajectories through a small Python function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd062a2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_g(traj,dr = 0.05):\n",
    "    '''\n",
    "    Calculates the 2D radial distribution function (RDF)\n",
    "    g(r) of a trajectory (list of ASE Atoms objects)\n",
    "    '''\n",
    "    Lx = traj[0].cell[0,0]\n",
    "    Ly = traj[0].cell[1,1]\n",
    "    #We create a list of Nx2 arrays\n",
    "    traj_2D = [a.positions[:,:2] for a in traj]\n",
    "    N = traj[0].get_global_number_of_atoms()\n",
    "    print(Lx,Ly,N)\n",
    "    def distance(dx, dy):\n",
    "        dx = dx - Lx*int(2/Lx*dx)\n",
    "        dy = dy - Ly*int(2/Ly*dy)\n",
    "        return dx**2 + dy**2\n",
    "\n",
    "    nmcs = len(traj_2D)\n",
    "    ngcum = 0\n",
    "    nbin = int(3.5/dr)\n",
    "    print('{} bins'.format(nbin))\n",
    "    gcum = np.zeros(nbin)\n",
    "\n",
    "    for sn in traj_2D:\n",
    "        x = sn[:,0] \n",
    "        y = sn[:,1]\n",
    "        for i in range(N-1):\n",
    "            for j in range(i+1,N):\n",
    "                dx = x[i] - x[j]\n",
    "                dy = y[i] - y[j]\n",
    "                r2 =  distance(dx,dy)\n",
    "                ibin = int(np.sqrt(r2)/dr) +1\n",
    "                if (ibin<nbin):\n",
    "                    gcum[ibin] = gcum[ibin] + 1 \n",
    "    rho = float(N)/(Lx*Ly)\n",
    "    xnorm =2.0/(rho*N*nmcs)\n",
    "    rl,gl = [],[]    \n",
    "    for ir in range(nbin):\n",
    "        r    = ir*dr + 0.5*dr  \n",
    "        area = 2.0*np.pi*r*dr    \n",
    "        g    = gcum[ir]*xnorm/area\n",
    "        rl.append(r)\n",
    "        gl.append(g)\n",
    "    return rl, gl\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee0e5298",
   "metadata": {},
   "outputs": [],
   "source": [
    "rl,gl = calculate_g(traj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "033a3778",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plotting\n",
    "plt.plot(rl,gl,label='Post-processing RDF')\n",
    "plt.xlabel('r')\n",
    "plt.ylabel('g(r)')\n",
    "#plt.plot(data[:,0],data[:,1])\n",
    "plt.xlim([0,3.3])\n",
    "plt.legend()\n",
    "plt.axhline(1,color='grey',alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fb2fecc",
   "metadata": {},
   "source": [
    "> Let's comparare it with the one calculated during the run by the Fortran code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f237c825",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.loadtxt('./Fortran_code/g_of_r.dat')\n",
    "#Plotting\n",
    "plt.plot(rl,gl,'-',label='Post-processing RDF')\n",
    "plt.xlabel('r')\n",
    "plt.ylabel('g(r)')\n",
    "plt.xlim([0,3.3])\n",
    "plt.plot(data[:,0],data[:,1], 'o',alpha=0.5,label='On-the-fly RDF')\n",
    "plt.legend()\n",
    "plt.axhline(1,color='grey',alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcd374e8",
   "metadata": {},
   "source": [
    "> Let's run some MC with the Fortran code `hd-MC_v2.x`for higher densities \n",
    "(N=16, 4.2x3.63, 500 steps for equilibration, 5000 for sampling)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74097b5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "traj2 = read('./Fortran_code/traj2.xyz',index=':')\n",
    "rl2,gl2 = calculate_g(traj2,dr=0.005)\n",
    "#data = np.loadtxt('./Fortran_code/g_of_r.dat')\n",
    "plt.plot(rl2,gl2,'-',label='Post-processing RDF')\n",
    "#plt.plot(data[:,0],data[:,1], 'o',alpha=0.5,label='On-the-fly RDF')\n",
    "plt.xlabel('r')\n",
    "plt.ylabel('g(r)')\n",
    "plt.xlim([0,3.3])\n",
    "plt.legend()\n",
    "plt.axhline(1,color='grey',alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f07e5595",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "sn = traj2[500]\n",
    "plot_atoms(sn,ax,radii=[0.5]*sn.get_global_number_of_atoms())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11eb59cf",
   "metadata": {},
   "source": [
    "> Let's run some MC with the Fortran code `hd-MC_v2.x`for intermediate densities \n",
    "(N=16, 4.45x3.85, 500 steps for equilibration, 5000 for sampling)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59272c63",
   "metadata": {},
   "outputs": [],
   "source": [
    "traj3 = read('./Fortran_code/traj3.xyz',index=':')\n",
    "rl3,gl3 = calculate_g(traj3,dr=0.005)\n",
    "plt.plot(rl,gl,'-',label=r'Post-processing RDF ($\\rho=0.49$)')\n",
    "plt.plot(rl2,gl2,'-',label=r'Post-processing RDF ($\\rho=0.91$)')\n",
    "plt.plot(rl3,gl3,'-',label=r'Post-processing RDF ($\\rho=0.81$)')\n",
    "plt.xlabel('r')\n",
    "plt.ylabel('g(r)')\n",
    "plt.xlim([0,2.5])\n",
    "plt.legend()\n",
    "plt.axhline(1,color='grey',alpha=0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30e85207",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "sn = traj3[500]\n",
    "plot_atoms(sn,ax,radii=[0.5]*sn.get_global_number_of_atoms())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54b69900",
   "metadata": {},
   "outputs": [],
   "source": [
    "view(traj3)"
   ]
  }
 ],
 "metadata": {
  "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.11.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}