{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import glob\n",
    "import MDAnalysis as mda\n",
    "from MDAnalysis.analysis.rms import rmsd\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as colors\n",
    "import matplotlib.ticker as ticker\n",
    "from matplotlib import cm\n",
    "from ipywidgets import *\n",
    "from matplotlib.widgets import Slider"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rmsd(posM):\n",
    "    # root mean square deviation in a MD trajectory\n",
    "    # return a rmsd array in 3 coordinates\n",
    "    initpos=posM[0]\n",
    "    lentraj=len(posM)\n",
    "    rmsd=np.zeros((lentraj))\n",
    "    for a in range(lentraj):\n",
    "        #rmsd[a]=np.sqrt(np.sum(np.square(initpos-posM[a]), axis=0)/np.shape(posM)[1])\n",
    "        rmsd[a]=np.sqrt(np.sum(np.square(initpos-posM[a]))/np.shape(posM)[1])\n",
    "    return rmsd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def msd(posM):\n",
    "    # mean square deviation in a MD trajectory\n",
    "    # return a rmsd array in 3 coordinates\n",
    "    initpos=posM[0]\n",
    "    lentraj=len(posM)\n",
    "    rmsd=np.zeros((lentraj))\n",
    "    for a in range(lentraj):\n",
    "        #rmsd[a]=np.sqrt(np.sum(np.square(initpos-posM[a]), axis=0)/np.shape(posM)[1])\n",
    "        rmsd[a]=np.sum(np.square(initpos-posM[a]))/np.shape(posM)[1]\n",
    "    return rmsd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def msdperatom(posM):\n",
    "    # mean square deviation in a MD trajectory\n",
    "    # return a rmsd array in 3 coordinates\n",
    "    initpos=posM[0]\n",
    "    lentraj=len(posM)\n",
    "    totatoms=np.shape(posM)[1]\n",
    "    rmsd=np.zeros((lentraj,totatoms))\n",
    "    for a in range(lentraj):\n",
    "        rmsd[a]=np.sum(np.square(initpos-posM[a]), axis=1)\n",
    "        #rmsd[a]=np.sqrt(np.sum(np.square(initpos-posM[a]))/np.shape(posM)[1])\n",
    "    return rmsd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def indexatmsel(atmsel1,atmsel2):\n",
    "    # return the initial and final index position where\n",
    "    # the elements of atmsel2 are in atmsel1\n",
    "    # atmsel1 > atmsel2\n",
    "    posini=np.where(atmsel1.indices == atmsel2.indices[0])[0][0]\n",
    "    posfin=np.where(atmsel1.indices == atmsel2.indices[-1])[0][0]\n",
    "    return posini,posfin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pbccorrect(rmsd):\n",
    "    dim1,dim2=np.shape(rmsd)\n",
    "    result=np.zeros((dim1,dim2))\n",
    "    for a in range(dim1):\n",
    "        for b in range(dim2):\n",
    "            if rmsd[a,b] > 100:\n",
    "                result[a,b]=rmsd[a,b]-174.894\n",
    "            else:\n",
    "                result[a,b]=rmsd[a,b]\n",
    "    return result        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def autocorr(x):\n",
    "    n=len(x)\n",
    "    ext_size=2*n-1\n",
    "    fsize=2**np.ceil(np.log2(ext_size)).astype('int')\n",
    "    xp=x-np.mean(x)\n",
    "    var=np.var(x)\n",
    "    cf=np.fft.fft(xp,fsize)\n",
    "    sf=cf.conjugate()*cf\n",
    "    corr=np.fft.ifft(sf).real\n",
    "    corr=corr/var/n\n",
    "    return corr[:len(x)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fftauto(vM):\n",
    "    global timestep\n",
    "    if len(np.shape(vM)) == 3:\n",
    "        coor=np.shape(vM)[2]\n",
    "        n_atom=np.shape(vM)[1]\n",
    "        vacf2=np.zeros((coor,len(vM)))\n",
    "        for dim in range(coor):\n",
    "            vacf3=np.zeros((n_atom,len(vM)))\n",
    "            for e in range(n_atom):\n",
    "                vacf3[e]=autocorr(vM[:,e,dim])\n",
    "            vacf2[dim]=vacf3.sum(axis=0)\n",
    "        vacf1=vacf2.sum(axis=0)\n",
    "    elif len(np.shape(vM)) == 2:\n",
    "        n_atom=np.shape(vM)[1]\n",
    "        vacf2=np.zeros((n_atom,len(vM)))\n",
    "        for e in range(n_atom):\n",
    "            vacf2[e]=autocorr(vM[:,e])\n",
    "        vacf1=vacf2.sum(axis=0)\n",
    "    else:\n",
    "        vacf1=autocorr(vM)\n",
    "    fftc1=np.fft.fft(vacf1)\n",
    "    fftm1=np.absolute(fftc1)\n",
    "    freq=np.fft.fftfreq(len(vM),d=timestep)\n",
    "    return freq[:len(vM)//2],fftm1[:len(vM)//2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 387,
   "metadata": {},
   "outputs": [],
   "source": [
    "#data=np.load(\"../charmm/dynapull/004/rmsdall_10.1.1.98res.npy\")\n",
    "#rmsdselall=pbccorrect(data)\n",
    "rmsdselall=np.load(\"../charmm/dynapull/free/1k/001/rmsdall_25000fr_miniabnr.98res.npy\")#RMSD data\n",
    "#velonormall=np.load(\"../charmm/dynapull/free/005/velnormpermol_98mol_10.1.0.vel.npy\")\n",
    "#fftselall=np.load(\"../charmm/dynapull/free/005/fftpermol2_98mol_10.1.0.vel.npy\")\n",
    "#coord=mda.Universe(\"../charmm/dynapull/inpcharmm.psf\",\"../charmm/dynapull/inpcharmm.crd\")\n",
    "#tempst=np.loadtxt(\"../charmm/dynapull/free/005/temps.dat\")\n",
    "try:\n",
    "    stt=np.where(tempst[:,0] == 4.)[0][0]-1\n",
    "    temps=tempst[stt:]\n",
    "except IndexError:\n",
    "    temps=tempst"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 364,
   "metadata": {},
   "outputs": [],
   "source": [
    "lentraj=len(rmsdselall)\n",
    "#strsel=\"resid 2 or resid 6 or resid 10 or resid 14 or resid 18 or resid 22 or resid 26 and not (name H*)\"\n",
    "strsel=\"resname DAB and name N3\"\n",
    "atmsel=coord.select_atoms(strsel)\n",
    "totatoms=atmsel.n_atoms\n",
    "lresids=np.unique(np.array(atmsel.resids))\n",
    "initsel=coord.select_atoms(strsel)\n",
    "initpos=initsel.positions\n",
    "X=initpos[:,0]\n",
    "Y=initpos[:,1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Kinetic energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 365,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 115916 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "plt.title(\"velnorm edge\\n\",fontsize=20)\n",
    "plt.ylabel('$|v_{mol}| (\\AA/ps$)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Time (ps)',fontsize=30,labelpad=15)\n",
    "colors=['tab:blue','tab:orange','tab:green','white','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "k=0\n",
    "for p in range(4,29,4):\n",
    "    ax.plot(np.array(range(steps))/250,velonormall[:steps,p//2-1]+k*0.8, linestyle='-' , color=colors[k], linewidth=2, label='resid '+str(p))\n",
    "    k+=1\n",
    "p,k=16,3\n",
    "ax.plot(np.array(range(steps))/250,velonormall[:steps,p//2-1]*1+k*0.8, color='tab:red', linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "\n",
    "plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,100) # Set the x axis limits\n",
    "plt.ylim(0,7) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(1))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(0.5))\n",
    "#ax.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
    "#plt.grid('both')\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 366,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 132492 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "plt.title(\"velnorm to Bulk_x\\n\",fontsize=20)\n",
    "plt.ylabel('$|v_{mol}| (\\AA/ps$)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Time (ps)',fontsize=30,labelpad=15)\n",
    "colors=['white','tab:blue','tab:orange','tab:green','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "k=0\n",
    "for p in range(16,185,28):\n",
    "    ax.plot(np.array(range(steps))/250,velonormall[:steps,p//2-1]+k*0.8, color=colors[k], linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "    k+=1\n",
    "p,k=16,0\n",
    "ax.plot(np.array(range(steps))/250,velonormall[:steps,p//2-1]*1+k*0.8, color='tab:red', linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "\n",
    "plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,100) # Set the x axis limits\n",
    "plt.ylim(0,7) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "#ax.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 367,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 115256 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "plt.title(\"velnorm other edge\\n\",fontsize=20)\n",
    "plt.ylabel('$|v_{mol}| (\\AA/ps$)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Time (ps)',fontsize=30,labelpad=15)\n",
    "colors=['tab:blue','tab:orange','tab:green','tab:olive','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "k=0\n",
    "for p in range(170,195,4):\n",
    "    ax.plot(np.array(range(steps))/250,velonormall[:steps,p//2-1]+k*0.8, color=colors[k], linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "    k+=1\n",
    "    #ax.plot(np.array(range(steps))/250,rmsdselall[:steps,12], color='purple', linestyle='-' , linewidth=5, label='resid 26')\n",
    "\n",
    "plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,100) # Set the x axis limits\n",
    "plt.ylim(0,7) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "#ax.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 368,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 129244 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "plt.title(\"velnorm 2nd-to-last edge\\n\",fontsize=20)\n",
    "plt.ylabel('$|v_{mol}| (\\AA/ps$)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Time (ps)',fontsize=30,labelpad=15)\n",
    "colors=['tab:blue','tab:orange','tab:green','tab:olive','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "k=0\n",
    "for p in range(2,27,4):\n",
    "    ax.plot(np.array(range(steps))/250,velonormall[:steps,p//2-1]+k*0.8, linestyle='-' , color=colors[k], linewidth=2, label='resid '+str(p))\n",
    "    k+=1\n",
    "#p,k=14,3\n",
    "#ax.plot(np.array(range(steps))/250,velonormall[:steps,p//2-1]*1+k*0.3, color='tab:red', linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "\n",
    "plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,100) # Set the x axis limits\n",
    "plt.ylim(0,7) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(1))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(0.5))\n",
    "#ax.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
    "#plt.grid('both')\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RMSD PLOTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 376,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 116252 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "plt.title(\"RMSD edge\\n\",fontsize=20)\n",
    "plt.ylabel('RMSD ($\\AA$)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Time (ps)',fontsize=30,labelpad=15)\n",
    "colors=['tab:blue','tab:orange','tab:green','white','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "k=0\n",
    "for p in range(4,29,4):\n",
    "    ax.plot(np.array(range(steps))/250,rmsdselall[:steps,p//2-1]+k*0.3, color=colors[k], linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "    k+=1\n",
    "p,k=16,3\n",
    "ax.plot(np.array(range(steps))/250,rmsdselall[:steps,p//2-1]*0.3+k*0.3, color='tab:red', linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "    \n",
    "    #ax.plot(np.array(range(steps))/250,rmsdselall[:steps,12], color='purple', linestyle='-' , linewidth=5, label='resid 26')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,10], color='blue', linestyle='-' , linewidth=5, label='resid 22')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,8], color='red', linestyle='-' , linewidth=5, label='resid 18')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,6], color='black', linestyle='-' , linewidth=5, label='resid 14')\n",
    "#ax.plot(rdf4[:,0], rdf4[:,1], color='purple', linestyle='-' , linewidth=5, label='rdf no surf')\n",
    "plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,100) # Set the x axis limits\n",
    "plt.ylim(0,2.5) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "#ax.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 375,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 126384 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "plt.title(\"RMSD to Bulk_x\\n\",fontsize=20)\n",
    "plt.ylabel('RMSD ($\\AA$)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Time (ps)',fontsize=30,labelpad=15)\n",
    "colors=['white','tab:blue','tab:orange','tab:green','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "k=0\n",
    "for p in range(16,185,28):\n",
    "    ax.plot(np.array(range(steps))/250,rmsdselall[:steps,p//2-1]+k*0.3, color=colors[k], linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "    k+=1\n",
    "p,k=16,0\n",
    "ax.plot(np.array(range(steps))/250,rmsdselall[:steps,p//2-1]*0.3+k*0.3, color='tab:red', linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "    #ax.plot(np.array(range(steps))/250,rmsdselall[:steps,12], color='purple', linestyle='-' , linewidth=5, label='resid 26')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,10], color='blue', linestyle='-' , linewidth=5, label='resid 22')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,8], color='red', linestyle='-' , linewidth=5, label='resid 18')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,6], color='black', linestyle='-' , linewidth=5, label='resid 14')\n",
    "#ax.plot(rdf4[:,0], rdf4[:,1], color='purple', linestyle='-' , linewidth=5, label='rdf no surf')\n",
    "plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,100) # Set the x axis limits\n",
    "plt.ylim(0,2.5) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "#ax.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 377,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 126988 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "plt.title(\"RMSD other edge\\n\",fontsize=20)\n",
    "plt.ylabel('RMSD ($\\AA$)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Time (ps)',fontsize=30,labelpad=15)\n",
    "colors=['tab:blue','tab:orange','tab:green','tab:olive','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "k=0\n",
    "for p in range(170,195,4):\n",
    "    ax.plot(np.array(range(steps))/250,rmsdselall[:steps,p//2-1]+k*0.3, color=colors[k], linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "    k+=1\n",
    "    #ax.plot(np.array(range(steps))/250,rmsdselall[:steps,12], color='purple', linestyle='-' , linewidth=5, label='resid 26')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,10], color='blue', linestyle='-' , linewidth=5, label='resid 22')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,8], color='red', linestyle='-' , linewidth=5, label='resid 18')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,6], color='black', linestyle='-' , linewidth=5, label='resid 14')\n",
    "#ax.plot(rdf4[:,0], rdf4[:,1], color='purple', linestyle='-' , linewidth=5, label='rdf no surf')\n",
    "plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,100) # Set the x axis limits\n",
    "plt.ylim(0,2.5) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "#ax.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 378,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 125068 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "plt.title(\"RMSD 2nd-to-last edge\\n\",fontsize=20)\n",
    "plt.ylabel('RMSD ($\\AA$)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Time (ps)',fontsize=30,labelpad=15)\n",
    "colors=['tab:blue','tab:orange','tab:green','tab:olive','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "k=0\n",
    "for p in range(2,27,4):\n",
    "    ax.plot(np.array(range(steps))/250,rmsdselall[:steps,p//2-1]+k*0.3, color=colors[k], linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "    k+=1\n",
    "    #ax.plot(np.array(range(steps))/250,rmsdselall[:steps,12], color='purple', linestyle='-' , linewidth=5, label='resid 26')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,10], color='blue', linestyle='-' , linewidth=5, label='resid 22')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,8], color='red', linestyle='-' , linewidth=5, label='resid 18')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,6], color='black', linestyle='-' , linewidth=5, label='resid 14')\n",
    "#ax.plot(rdf4[:,0], rdf4[:,1], color='purple', linestyle='-' , linewidth=5, label='rdf no surf')\n",
    "plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,100) # Set the x axis limits\n",
    "plt.ylim(0,2.5) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "#ax.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FFT of RMSD per molecule"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### FFT PLOTS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 379,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 90220 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "plt.title(\"FFT edge\\n\",fontsize=20)\n",
    "plt.ylabel('FFT ($10^{-3}$)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Frequency (THz)',fontsize=30,labelpad=15)\n",
    "colors=['tab:blue','tab:orange','tab:green','tab:red','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "k=0\n",
    "for p in range(4,29,4):\n",
    "    ax.plot(np.real(fftselall[0]),np.real(fftselall[p//2,:])/1000+k*0.6, color=colors[k], linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "    k+=1\n",
    "    #ax.plot(np.array(range(steps))/250,rmsdselall[:steps,12], color='purple', linestyle='-' , linewidth=5, label='resid 26')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,10], color='blue', linestyle='-' , linewidth=5, label='resid 22')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,8], color='red', linestyle='-' , linewidth=5, label='resid 18')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,6], color='black', linestyle='-' , linewidth=5, label='resid 14')\n",
    "#ax.plot(rdf4[:,0], rdf4[:,1], color='purple', linestyle='-' , linewidth=5, label='rdf no surf')\n",
    "plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,100) # Set the x axis limits\n",
    "plt.ylim(0,5) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "#ax.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 380,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 68748 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "plt.title(\"FFT to Bulk_x\\n\",fontsize=20)\n",
    "plt.ylabel('FFT ($10^{-3}$)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Frequency (THz)',fontsize=30,labelpad=15)\n",
    "colors=['tab:red','tab:blue','tab:orange','tab:green','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "k=0\n",
    "for p in range(16,185,28):\n",
    "    ax.plot(np.real(fftselall[0]),np.real(fftselall[p//2,:])/1000+k*1, color=colors[k], linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "    k+=1\n",
    "    #ax.plot(np.array(range(steps))/250,rmsdselall[:steps,12], color='purple', linestyle='-' , linewidth=5, label='resid 26')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,10], color='blue', linestyle='-' , linewidth=5, label='resid 22')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,8], color='red', linestyle='-' , linewidth=5, label='resid 18')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,6], color='black', linestyle='-' , linewidth=5, label='resid 14')\n",
    "#ax.plot(rdf4[:,0], rdf4[:,1], color='purple', linestyle='-' , linewidth=5, label='rdf no surf')\n",
    "#plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,100) # Set the x axis limits\n",
    "plt.ylim(0,7) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "#ax.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 381,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 73872 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "plt.title(\"FFT other edge\\n\",fontsize=20)\n",
    "plt.ylabel('FFT ($10^{-3}$)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Frequency (THz)',fontsize=30,labelpad=15)\n",
    "colors=['tab:blue','tab:orange','tab:green','tab:olive','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "k=0\n",
    "for p in range(170,195,4):\n",
    "    ax.plot(np.real(fftselall[0]),np.real(fftselall[p//2,:])/1000+k*1, color=colors[k], linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "    k+=1\n",
    "    #ax.plot(np.array(range(steps))/250,rmsdselall[:steps,12], color='purple', linestyle='-' , linewidth=5, label='resid 26')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,10], color='blue', linestyle='-' , linewidth=5, label='resid 22')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,8], color='red', linestyle='-' , linewidth=5, label='resid 18')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,6], color='black', linestyle='-' , linewidth=5, label='resid 14')\n",
    "#ax.plot(rdf4[:,0], rdf4[:,1], color='purple', linestyle='-' , linewidth=5, label='rdf no surf')\n",
    "#plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,100) # Set the x axis limits\n",
    "#plt.ylim(0,2.65) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "#ax.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 382,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 69472 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "plt.title(\"FFT 2nd-to-last edge\\n\",fontsize=20)\n",
    "plt.ylabel('FFT ($10^{-3}$)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Frequency (THz)',fontsize=30,labelpad=15)\n",
    "colors=['tab:blue','tab:orange','tab:green','tab:olive','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "k=0\n",
    "for p in range(2,27,4):\n",
    "    ax.plot(np.real(fftselall[0]),np.real(fftselall[p//2,:])/1000+k*1, color=colors[k], linestyle='-' , linewidth=2, label='resid '+str(p))\n",
    "    k+=1\n",
    "    #ax.plot(np.array(range(steps))/250,rmsdselall[:steps,12], color='purple', linestyle='-' , linewidth=5, label='resid 26')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,10], color='blue', linestyle='-' , linewidth=5, label='resid 22')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,8], color='red', linestyle='-' , linewidth=5, label='resid 18')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,6], color='black', linestyle='-' , linewidth=5, label='resid 14')\n",
    "#ax.plot(rdf4[:,0], rdf4[:,1], color='purple', linestyle='-' , linewidth=5, label='rdf no surf')\n",
    "#plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,100) # Set the x axis limits\n",
    "#plt.ylim(0,2.65) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "#ax.yaxis.set_major_locator(ticker.MultipleLocator(5))\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 383,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 24764 characters]\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the data\n",
    "steps=25000\n",
    "fig, ax = plt.subplots(figsize=(8,5)) # Figure size\n",
    "plt.title(\"Temperature\",fontsize=20)\n",
    "plt.ylabel('Temperature (K)',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Time (ps)',fontsize=30,labelpad=15)\n",
    "colors=['tab:red','tab:blue','tab:orange','tab:green','tab:purple',\n",
    "        'tab:brown','tab:pink','tab:gray','tab:olive','tab:cyan']\n",
    "\n",
    "ax.plot(temps[:,0]/1000,temps[:,1], color='tab:red', linestyle='-' , linewidth=2)#, label='resid '+str(p))\n",
    "\n",
    "    #ax.plot(np.array(range(steps))/250,rmsdselall[:steps,12], color='purple', linestyle='-' , linewidth=5, label='resid 26')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,10], color='blue', linestyle='-' , linewidth=5, label='resid 22')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,8], color='red', linestyle='-' , linewidth=5, label='resid 18')\n",
    "#ax.plot(np.array(range(steps))/250,rmsdselall[:steps,6], color='black', linestyle='-' , linewidth=5, label='resid 14')\n",
    "#ax.plot(rdf4[:,0], rdf4[:,1], color='purple', linestyle='-' , linewidth=5, label='rdf no surf')\n",
    "#plt.legend(fontsize=10)\n",
    "ax.set_xscale(\"log\")\n",
    "plt.xlim(0.1,10) # Set the x axis limits\n",
    "plt.ylim(0,20) # Set the y axis limits\n",
    "plt.xticks(fontsize=30)\n",
    "plt.yticks(fontsize=30)\n",
    "#ax.xaxis.set_major_formatter(ticker.FormatStrFormatter('%0.1f'))\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(50))\n",
    "#ax.xaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "ax.yaxis.set_major_locator(ticker.MultipleLocator(10))\n",
    "for axis in ['top','bottom','left','right']:\n",
    "    ax.spines[axis].set_linewidth(2.0)\n",
    "ax.xaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "ax.xaxis.set_tick_params(width=1.5,length=5,which='minor')\n",
    "ax.yaxis.set_tick_params(width=2.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.savefig(\"potential.png\") # The graphs can be as PNG, SVG, BMP or PDF"
   ]
  },
  {
   "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 (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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
