{
 "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 import lineardensity as lin\n",
    "#import nglview as nv\n",
    "#import ipywidgets\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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dcdtovel(u,sel,nm=1):\n",
    "    \"\"\" u and sel are MDAnalysis\n",
    "    universe and atomic selection.\n",
    "    returns the velocity vector\n",
    "    of the nm mode for the selection\"\"\"\n",
    "    posit=np.zeros((12, len(sel), 3))\n",
    "    k=0\n",
    "    for ts in u.trajectory[(nm-1)*12:nm*12]:\n",
    "        posit[k]=sel.positions\n",
    "        k+=1\n",
    "    velo=np.zeros((12,len(sel),3))\n",
    "    for at in range(len(sel)):        \n",
    "        velo[:,at]=np.gradient(posit[:,at], axis=0)\n",
    "    return velo[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dcd2omega(u,sel,nm=1):\n",
    "    \"\"\" u and sel are MDAnalysis\n",
    "    universe and atomic selection.\n",
    "    returns the velocity vector\n",
    "    of the nm mode for the selection\"\"\"\n",
    "    enedab=u.select_atoms(\"type 9 and resname DAB\")\n",
    "    posMres=np.zeros((12,len(sel),3))\n",
    "    enevecs=np.zeros((enedab.n_residues,3))\n",
    "    posit=np.zeros((12, len(sel), 3))\n",
    "    atmperres=np.shape(posit)[1]//enedab.n_residues\n",
    "    k=0\n",
    "    for ts in u.trajectory[(nm-1)*12:nm*12]:\n",
    "        posit[k]=sel.positions\n",
    "        k+=1\n",
    "    for e in range(enedab.n_residues):\n",
    "        enevecs[e]=sel[atmperres*e:atmperres*(e+1)].center_of_mass()\n",
    "        posMres[:,atmperres*e:atmperres*(e+1)]=posit[:,atmperres*e:atmperres*(e+1)]-enevecs[e]\n",
    "    velo=np.zeros((12,len(sel),3))\n",
    "    for at in range(len(sel)):        \n",
    "        velo[:,at]=np.gradient(posit[:,at], axis=0)\n",
    "    rsqr=np.linalg.norm(posMres,axis=2)**(-2)\n",
    "    invpos=np.zeros((np.shape(posMres)))\n",
    "    for trj in range(np.shape(posMres)[0]):\n",
    "        for at in range(np.shape(posMres)[1]):\n",
    "            invpos[trj,at]=posMres[trj,at]*rsqr[trj,at]\n",
    "    ome=np.cross(invpos,velo)\n",
    "    return ome[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mergeeigfrq(entry_list,dec=1):\n",
    "    # groups the the eigenvalues/NM within a certain deltaFreq\n",
    "    # return a list of list with initial and final limits for each mode\n",
    "    incr=np.power(10.,dec)\n",
    "    round_list=np.around(entry_list*incr)\n",
    "    nestM=[]\n",
    "    count=0\n",
    "    while count < len(round_list):\n",
    "        nest1=[]\n",
    "        nest1.append(count+1)\n",
    "        ref_value=round_list[count]\n",
    "        try:\n",
    "            next_value=round_list[count+1] # next value\n",
    "        except IndexError: # account for last value\n",
    "            nest1.append(count+1)\n",
    "            nestM.append(nest1)\n",
    "            break\n",
    "        while ref_value == next_value:\n",
    "            count+=1\n",
    "            try:\n",
    "                next_value=round_list[count+1] ### problem!!!!\n",
    "            except IndexError:\n",
    "                break#nextf+=deltf\n",
    "        nest1.append(count+1) # final index\n",
    "        nestM.append(nest1)\n",
    "        count+=1\n",
    "    return nestM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def intensnm(NM,minf):\n",
    "    # return a list with the NM with intensities above minf\n",
    "    intensnm=NM[~np.any(NM <= 0, axis=1)].T\n",
    "    intensnm1=np.vstack((np.arange(1,len(intensnm.T)+1),intensnm[0],intensnm[1]))\n",
    "    intensnm2=intensnm1.T[~np.any(intensnm1.T <= minf, axis=1)].T\n",
    "    return intensnm2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mergebigfrq(entry_list,NM,minf,dec=1):\n",
    "    intpeaks=np.asarray(intensnm(NMspec,minf)[0]).astype(int)\n",
    "    bigfreq=[]\n",
    "    for p in range(len(intpeaks)):\n",
    "        ix=intpeaks[p]\n",
    "        bigfreq.append(eigfreq[ix])\n",
    "    bigfreq=np.asarray(bigfreq)\n",
    "    frecum=mergeeigfrq(bigfreq,dec)\n",
    "    return intpeaks,frecum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Input data for mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load NM trajectories, coordinates, NMA spectra and eigenfreq\n",
    "# Ladder SSH Topological free boundary case\n",
    "nmtrajfree=mda.Universe(\"charmm/inpc.copfe.52.psf\",\"charmm/z-ingap01/nmodes.copfe.52.free.dcd\")\n",
    "coord=mda.Universe(\"charmm/inpc.copfe.52.psf\",\"charmm/z-ingap01/minisd.copfe.52.pdb\")\n",
    "NMspecfree=np.loadtxt(\"charmm/z-ingap01/NMAspectra_nma.52.free.COP.dat\")\n",
    "eigfreqfree=np.loadtxt(\"charmm/z-ingap01/eigenfrq_nma.52.free.COP.dat\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load NM trajectories, coordinates, NMA spectra and eigenfreq\n",
    "# Ladder SSH Topological free heavy boundary case\n",
    "nmtrajfreem=mda.Universe(\"charmm/inpc.copfe.52.psf\",\"charmm/z-ingap01/nmodes.copfe.52.freem.dcd\")\n",
    "coord=mda.Universe(\"charmm/inpc.copfe.52.psf\",\"charmm/z-ingap01/minisd.copfe.52.pdb\")\n",
    "NMspecfreem=np.loadtxt(\"charmm/z-ingap01/NMAspectra_nma.52.freem.COP.dat\")\n",
    "eigfreqfreem=np.loadtxt(\"charmm/z-ingap01/eigenfrq_nma.52.freem.COP.dat\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load NM trajectories, coordinates, NMA spectra and eigenfreq\n",
    "# Ladder SSH Topological free boundary strong boundary case\n",
    "nmtrajfreestr=mda.Universe(\"charmm/inpc.copfe.52.psf\",\"charmm/z-ingap01/nmodes.copfe.52.freestr.dcd\")\n",
    "coord=mda.Universe(\"charmm/inpc.copfe.52.psf\",\"charmm/z-ingap01/minisd.copfe.52.pdb\")\n",
    "NMspecfreestr=np.loadtxt(\"charmm/z-ingap01/NMAspectra_nma.52.freestr.COP.dat\")\n",
    "eigfreqfreestr=np.loadtxt(\"charmm/z-ingap01/eigenfrq_nma.52.freestr.COP.dat\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "#nmtraj,eigfreq,NMspec=nmtrajfree,eigfreqfree,NMspecfree\n",
    "nmtraj,eigfreq,NMspec=nmtrajfreem,eigfreqfreem,NMspecfreem\n",
    "#nmtraj,eigfreq,NMspec=nmtrajfreestr,eigfreqfreestr,NMspecfreestr\n",
    "#nmtraj,eigfreq,NMspec=nmtrajpbc,eigfreqpbc,NMspecpbc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "#mergefrecu=mergeeigfrq(eigfreq,-0.1)\n",
    "bigfreq,mergefrecu=mergebigfrq(eigfreq,NMspec,0.05,0.5)\n",
    "lentraj=len(nmtraj.trajectory)//12\n",
    "leneigfrq=len(eigfreq)\n",
    "#assert (lentraj == leneigfrq), \"Trajectory movie and eigenfrequency count are not the same length\"\n",
    "strsel=\"bynum 1:52\" # atomic system selection for the plot\n",
    "atmsel=nmtraj.select_atoms(strsel)\n",
    "initsel=coord.select_atoms(strsel)\n",
    "initpos=initsel.positions\n",
    "X=initpos[:,0]\n",
    "Y=initpos[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 18668 characters]\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### Plot the spectrum/spectra\n",
    "fig, ax = plt.subplots(figsize=(8,8)) # Figure size\n",
    "#plt.title(\"NMA crystal and Nanoribbon\")\n",
    "#plt.ylabel('Intensity',fontsize=30,labelpad=10)\n",
    "plt.xlabel('Frequency (cm$^{-1}$)',fontsize=20,labelpad=15)\n",
    "\n",
    "#ax.plot(NMspec[:,0], NMspec[:,1], color='tab:gren' , label='Nanoribbon')\n",
    "ax.plot(NMspec[:,0], NMspec[:,1], color='black' )#, label='Hex 73') \n",
    "#ax.scatter(NMspec[:,0], NMspec[:,1], color='black' )#, label='Hex 73') \n",
    "\n",
    "#ax.plot(NMspec2[:,0], NMspec2[:,1]*-1+2.5, color='tab:red' , label='Crystal')\n",
    "\n",
    "#plt.xlim(-5,600) #Set the x axis limits\n",
    "#plt.xscale(\"log\")\n",
    "plt.ylim(0,2.5) # Set the y axis limits\n",
    "plt.xticks(fontsize=20)\n",
    "plt.yticks(fontsize=10)\n",
    "#ax.yaxis.set_major_formatter(ticker.FormatStrFormatter(''))\n",
    "ax.xaxis.set_major_locator(ticker.MultipleLocator(100))\n",
    "ax.xaxis.set_minor_locator(ticker.MultipleLocator(25))\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=0.0,length=10,which='major')\n",
    "plt.tight_layout()\n",
    "#plt.legend(fontsize=15,loc=0)\n",
    "plt.grid(which='both')\n",
    "#plt.savefig(\"NMA.pdp.112.aa.80-110.png\") # The last spectra can be saved as an image file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quiver or scatter plot based on initial and final normal mode velocities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 19616 characters]\n",
      "text/plain": [
       "<Figure size 1440x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#mm=2\n",
    "#fini,ffin=bigfreq[mergefrecu[mm-1][:2][0]-1],bigfreq[mergefrecu[mm-1][:2][1]-1]\n",
    "fck=288\n",
    "fini,ffin=fck,fck # initial normal mode, final normal mode\n",
    "#plt.figure(figsize=(20,4))#,(8.75,7.5)\n",
    "fig, ax = plt.subplots(figsize=(20,5))\n",
    "#plt.figure(figsize=(4,20))\n",
    "ff=np.average(eigfreq[fini-1:ffin])\n",
    "#plt.title('ROC Curve; Enrichment factor 1% = {:.2f}'.format(ff))\n",
    "#plt.title('eigfrequency: {:.2f} cm-1\\n Modes: {:d} - {:d}'.format(ff,fini,ffin))\n",
    "velonm=np.zeros((ffin+1-fini,len(atmsel),3))\n",
    "for nm in range(fini,ffin+1):\n",
    "    velonm[nm-fini]=dcdtovel(nmtraj,atmsel,nm)\n",
    "velonma=np.sum(velonm,axis=0)\n",
    "veloabs=np.linalg.norm(velonma,axis=1)\n",
    "velabsn=veloabs/np.power(np.max(veloabs),2)\n",
    "veloz=np.abs(velonma[:,2])\n",
    "if np.abs(np.min(velonma[:,0])/np.max(velonma[:,1])) > 10000:\n",
    "    velx,vely=velonma[:,1],velonma[:,0]\n",
    "    plt.text(0, -45, 'Longitudinal', fontsize = 30)\n",
    "else:\n",
    "    velx,vely=velonma[:,0],velonma[:,1]\n",
    "    plt.text(-10, -45, 'Tranversal', fontsize = 23)\n",
    "ax.quiver(X,Y, velx,vely, scale_units='xy', scale=np.max(velonma[:,:2])/55, cmap=\"jet\",linewidth=1)\n",
    "#    ax=plt.scatter(X,Y, c=veloabs, s=velabsn, cmap=\"rainbow\")\n",
    "#plt.xlim(-90,100)\n",
    "#plt.xlim(-60,80)\n",
    "plt.xlim(-10,110)\n",
    "#plt.ylim(-65,85)\n",
    "plt.ylim(-60,60)\n",
    "plt.xticks(fontsize=30)\n",
    "#plt.yticks(fontsize=0)\n",
    "ax.yaxis.set_ticks([])\n",
    "cb.ax.tick_params(labelsize=15)\n",
    "#plt.savefig(\"nmode\"+str(fini)+\"-\"+str(ffin)+\".pdbr.96.nr.ave.png\",transparent=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.15480877459049225, -207646125327413.56)"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(velonma[:,:3]),np.min(velonma[:,0])/np.max(velonma[:,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(52, 3)"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.shape(velonma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/4m/nrl0c6_144v6spljvv1x4_y00000gq/T/ipykernel_27632/694029008.py:3: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  plt.figure(figsize=(20,4))#,(8.75,7.5)\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 39172 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 39988 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 39384 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 38916 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 37360 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 38776 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34408 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 38944 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 33476 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 32540 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 40820 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 43452 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 41080 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 41392 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 41940 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 41784 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 43012 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 43640 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 42932 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 44076 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 45004 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 43588 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 45392 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 41396 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 41008 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 40980 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 43036 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 40884 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 40796 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 27640 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 28816 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 29920 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 30264 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 36556 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 28516 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 28532 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 37092 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 29964 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 31696 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 35312 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 30712 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 37696 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 34108 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 30560 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 37536 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 36212 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 31888 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 32260 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 32084 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 32440 characters]\n",
      "text/plain": [
       "<Figure size 1440x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for i in range(200,250):\n",
    "    fini,ffin=i,i # initial normal mode, final normal mode\n",
    "    plt.figure(figsize=(20,4))#,(8.75,7.5)\n",
    "    #plt.figure(figsize=(4,20))\n",
    "    ff=np.average(eigfreq[fini-1:ffin])\n",
    "    #plt.title('ROC Curve; Enrichment factor 1% = {:.2f}'.format(ff))\n",
    "    plt.title('eigfrequency: {:.2f} cm-1\\n Modes: {:d} - {:d}'.format(ff,fini,ffin))\n",
    "    velonm=np.zeros((ffin+1-fini,len(atmsel),3))\n",
    "    for nm in range(fini,ffin+1):\n",
    "        velonm[nm-fini]=dcdtovel(nmtraj,atmsel,nm)\n",
    "    velonma=np.average(velonm,axis=0)\n",
    "    veloabs=np.linalg.norm(velonma,axis=1)\n",
    "    velabsn=veloabs/np.power(np.max(veloabs),2)\n",
    "    veloz=np.abs(velonma[:,2])\n",
    "    if np.abs(np.min(velonma[:,0])/np.max(velonma[:,1])) > 10000:\n",
    "        velx,vely=velonma[:,1],velonma[:,0]\n",
    "        plt.text(-10, -15, 'Longitudinal', fontsize = 23)\n",
    "    else:\n",
    "        velx,vely=velonma[:,0],velonma[:,1]\n",
    "        plt.text(-10, -15, 'Tranversal', fontsize = 23)\n",
    "    ax=plt.quiver(X,Y, velx,vely, veloz, scale_units='xy', scale=np.max(velonma[:,:2])/10, cmap=\"jet\",linewidth=1)\n",
    "    #    ax=plt.scatter(X,Y, c=veloabs, s=velabsn, cmap=\"rainbow\")\n",
    "    #plt.xlim(-90,100)\n",
    "    #plt.xlim(-60,80)\n",
    "    plt.xlim(-30,130)\n",
    "    #plt.ylim(-65,85)\n",
    "    plt.ylim(-20,20)\n",
    "    cb = plt.colorbar(ax)\n",
    "    cb.ax.tick_params(labelsize=15)\n",
    "    #plt.savefig(\"nmode\"+str(fini)+\"-\"+str(ffin)+\".pdbr.96.nr.ave.png\",transparent=True)"
   ]
  },
  {
   "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
}
