{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0c63f72a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "import numpy as np\n",
    "import glob\n",
    "import MDAnalysis as mda\n",
    "#from functools import partial\n",
    "#import pathos.multiprocessing as mp\n",
    "#import time\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": 56,
   "id": "b7a0cc61",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_nearest(array, value):\n",
    "    array = np.asarray(array)\n",
    "    idx = (np.abs(array - value)).argmin()\n",
    "    return idx"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef569483",
   "metadata": {},
   "source": [
    "### Copfe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "id": "cd365644",
   "metadata": {},
   "outputs": [],
   "source": [
    "## In-bulk topological kappa_1,_2= 3000,1000 C-C-C-C-C... kappa_3 = kappa_4 = 1000\n",
    "## folder copfe/008, similar to copfe/004 but with b_0 = 2Å for all bonds\n",
    "## Figure 2 paper!!\n",
    "psdfreeK=np.load(\"./gromacs/copfe/008/psdperatom_nvt.copfe.52.freeK.nj.52at.npy\") # k_strong=10^6kJ/mol/nm\n",
    "psdfreeM=np.load(\"./gromacs/copfe/008/psdperatom_nvt.copfe.52.freeM.nj.52at.npy\")\n",
    "psdpbc=np.load(\"./gromacs/copfe/008/psdperatom_nvt.copfe.52.pbc.nj.52at.npy\")\n",
    "psdfree=np.load(\"./gromacs/copfe/008/psdperatom_nvt.copfe.52.free.nj.52at.npy\")\n",
    "psdfix=np.load(\"./gromacs/copfe/008/psdperatom_nvt.copfe.52.fix.nj.52at.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "id": "e50e2616",
   "metadata": {},
   "outputs": [],
   "source": [
    "## In-bulk topological kappa_1,_2= 3000,1000 C-C-C-C-C... kappa_3 = kappa_4 = 1000\n",
    "## folder copfe/012, similar to copfe/004 but with b_0 = 0Å for all bonds\n",
    "psdfreeK=np.load(\"./gromacs/copfe/012/psdperatom_nvt.copfe.52.freeK.nj.52at.npy\") # k_strong=10^6kJ/mol/nm\n",
    "psdfreeM=np.load(\"./gromacs/copfe/012/psdperatom_nvt.copfe.52.freeM.nj.52at.npy\")\n",
    "psdpbc=np.load(\"./gromacs/copfe/012/psdperatom_nvt.copfe.52.pbc.nj.52at.npy\")\n",
    "psdfree=np.load(\"./gromacs/copfe/012/psdperatom_nvt.copfe.52.free.nj.52at.npy\")\n",
    "psdfix=np.load(\"./gromacs/copfe/012/psdperatom_nvt.copfe.52.fix.nj.52at.npy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11831639",
   "metadata": {},
   "source": [
    "### Ladder SSH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "id": "0fd883ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ladder SSH, topological case, kappa1,kappa2,kappa3=3000,1000,3000, zigzag 60° substrate, no restriction in main chain\n",
    "# folder ladderSSH/020\n",
    "# equil distance to substrate 0Å\n",
    "#psdfree=np.load(\"./gromacs/ladderSSH/020/psdperatom_nvt.copfe.52.free.nj.52at.npy\")\n",
    "#psdfreeK=np.load(\"./gromacs/ladderSSH/020/psdperatom_nvt.copfe.52.freeK.nj.52at.npy\")\n",
    "psdfreeM=np.load(\"./gromacs/ladderSSH/020/psdperatom_nvt.copfe.52.freeM.nj.52at.npy\")\n",
    "#psdpbc=np.load(\"./gromacs/ladderSSH/020/psdperatom_nvt.copfe.52.pbc.nj.52at.npy\")\n",
    "#psdfix=np.load(\"./gromacs/ladderSSH/020/psdperatom_nvt.copfe.52.fix.nj.52at.npy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71e927eb",
   "metadata": {},
   "source": [
    "### Triangular SSH ladder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "id": "1df183dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "#triangular laddder, topological case, kappa1,kappa2,kappa4,kappa5=3000,1000,500,500\n",
    "#ladder kSz(kappa3)=3000kJ/molnm\n",
    "#ladder dSz=0Å\n",
    "# folder triladder/003\n",
    "psdfree=np.load(\"./gromacs/triladder/003/psdperatom_nvt.trilad.52.free.nj.52at.npy\")\n",
    "psdfreeK=np.load(\"./gromacs/triladder/003/psdperatom_nvt.trilad.52.freeK.nj.52at.npy\")\n",
    "psdfreeM=np.load(\"./gromacs/triladder/003/psdperatom_nvt.trilad.52.freeM.nj.52at.npy\")\n",
    "psdpbc=np.load(\"./gromacs/triladder/003/psdperatom_nvt.trilad.52.pbc.nj.52at.npy\")\n",
    "psdfix=np.load(\"./gromacs/triladder/003/psdperatom_nvt.trilad.52.fix.nj.52at.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "id": "24588d28",
   "metadata": {},
   "outputs": [],
   "source": [
    "#triangular laddder, topological case, kappa1,kappa2,kappa4,kappa5=3000,2000,250,250\n",
    "#ladder kSz(kappa3)=1500kJ/molnm\n",
    "#ladder dSz=0Å\n",
    "# folder triladder/012\n",
    "psdfree=np.load(\"./gromacs/triladder/012/psdperatom_nvt.trilad.52.free.nj.52at.npy\")\n",
    "psdfreeK=np.load(\"./gromacs/triladder/012/psdperatom_nvt.trilad.52.freeK.nj.52at.npy\")\n",
    "psdfreeM=np.load(\"./gromacs/triladder/012/psdperatom_nvt.trilad.52.freeM.nj.52at.npy\")\n",
    "psdpbc=np.load(\"./gromacs/triladder/012/psdperatom_nvt.trilad.52.pbc.nj.52at.npy\")\n",
    "psdfix=np.load(\"./gromacs/triladder/012/psdperatom_nvt.trilad.52.fix.nj.52at.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "id": "9397ed44",
   "metadata": {},
   "outputs": [],
   "source": [
    "#triangular laddder, trivial case, kappa1,kappa2,kappa4,kappa5=2000,3000,250,250\n",
    "#ladder kSz(kappa3)=1500kJ/molnm\n",
    "#ladder dSz=0Å\n",
    "# folder triladder/012/trivial\n",
    "#psdfree=np.load(\"./gromacs/triladder/012/trivial/psdperatom_nvt.trilad.52.free.nj.52at.npy\")\n",
    "#psdfreeK=np.load(\"./gromacs/triladder/012/trivial/psdperatom_nvt.trilad.52.freeK.nj.52at.npy\")\n",
    "psdfreeM=np.load(\"./gromacs/triladder/012/trivial/psdperatom_nvt.trilad.52.freeM.nj.52at.npy\")\n",
    "#psdpbc=np.load(\"./gromacs/triladder/012/trivial/psdperatom_nvt.trilad.52.pbc.nj.52at.npy\")\n",
    "#psdfix=np.load(\"./gromacs/triladder/012/trivial/psdperatom_nvt.trilad.52.fix.nj.52at.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "daa05b5d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.5, 0. , 0. , 1. ])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "col"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 249,
   "id": "d1ef480a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 41040 characters]\n",
      "text/plain": [
       "<Figure size 576x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# PLOT of FFT[e] of different pairs of pearls\n",
    "fig,ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "dataplot=psdfreeM\n",
    "X=dataplot[0]/(33.35641)#*np.sqrt(0.4184))\n",
    "Y=dataplot[1:]\n",
    "lenY=len(Y)-19\n",
    "ixi,ixf=find_nearest(np.real(X), 0.28),find_nearest(np.real(X), 0.295)\n",
    "ymax=[]\n",
    "xymax=[]\n",
    "#plt.figure(figsize=(10,15))\n",
    "#plt.title(\"PSD away from boundary\")\n",
    "colors = iter(cm.jet(np.linspace(0, 1, lenY)))\n",
    "for e in range(0,lenY,1):\n",
    "    yvalue=0.01*np.real(Y[e])/np.average(np.real(Y))\n",
    "    ymax.append(np.max(yvalue[ixi:ixf])*0.5+7.5)\n",
    "    xymax.append(0.35+0.0027*e)\n",
    "    #col = next(colors)\n",
    "    col = ['black','blue','tab:cyan','green','yellow','orange','red']\n",
    "    ax.plot(np.real(X), yvalue-e+6,c=col[e])#, label=\"Ang = {:.1f}\".format(angles[a]))#,colors[e])\n",
    "#for e in range(len(Y)-15):\n",
    "#    plt.plot(np.real(X), 0.15*np.real(Y[e])/np.average(np.real(Y))+e,'b-')\n",
    "    #plt.plot(np.real(ffk[0][len(ffk[0])//2:]), np.real(ffk[e][len(ffk[0])//2:])/50000+e,'b-')\n",
    "#ax.plot(xymax,ymax, color='tab:red')#,linestyle=\"--\")\n",
    "plt.ylabel('PSD (Intensity)', fontsize=30)\n",
    "plt.xlabel('Frequency (cm$^{-1}$)', fontsize=30)\n",
    "plt.xlim(0.38,0.40)\n",
    "plt.ylim(0,11)\n",
    "plt.grid()\n",
    "plt.xticks(fontsize=22)\n",
    "plt.yticks(fontsize=22)\n",
    "ax.xaxis.set_major_locator(ticker.MultipleLocator(0.005))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "id": "81802062",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-0.5,  0. ,  0.5,  1. ,  1.5,  2. ,  2.5,  3. ]),\n",
       " [Text(0, 0, ''),\n",
       "  Text(0, 0, ''),\n",
       "  Text(0, 0, ''),\n",
       "  Text(0, 0, ''),\n",
       "  Text(0, 0, ''),\n",
       "  Text(0, 0, ''),\n",
       "  Text(0, 0, ''),\n",
       "  Text(0, 0, '')])"
      ]
     },
     "execution_count": 246,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 51884 characters]\n",
      "text/plain": [
       "<Figure size 612x720 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# PLOT of FFT[e] of different pairs of pearls\n",
    "fig,ax = plt.subplots(figsize=(8.5,10)) # Figure size\n",
    "dataplot=psdfreeM\n",
    "X=dataplot[0]/(33.35641)#*np.sqrt(0.4184))\n",
    "Y=dataplot[1:]\n",
    "lenY=len(Y)\n",
    "ixi,ixf=find_nearest(np.real(X), 0.37),find_nearest(np.real(X), 0.39)\n",
    "ymax=[]\n",
    "xymax=[]\n",
    "#plt.figure(figsize=(10,15))\n",
    "#plt.title(\"PSD away from boundary\")\n",
    "#colors = iter(cm.jet(np.linspace(0, 1, lenY)))\n",
    "for e in range(0,lenY,1):\n",
    "    yvalue=0.01*np.real(Y[e])/np.average(np.real(Y))\n",
    "    runy=np.convolve(yvalue, np.ones(10)/10, mode='same')\n",
    "    ymax.append(np.max(yvalue[ixi:ixf]))\n",
    "    #ymax.append(np.max(runy[ixi:ixf]))\n",
    "    xymax.append(e)\n",
    "    #col = next(colors)\n",
    "    #col = ['black','blue','tab:cyan','green','yellow','orange','red']\n",
    "    #ax.plot(np.real(X), yvalue-e+6,c=col[e])#, label=\"Ang = {:.1f}\".format(angles[a]))#,colors[e])\n",
    "#for e in range(len(Y)-15):\n",
    "#    plt.plot(np.real(X), 0.15*np.real(Y[e])/np.average(np.real(Y))+e,'b-')\n",
    "    #plt.plot(np.real(ffk[0][len(ffk[0])//2:]), np.real(ffk[e][len(ffk[0])//2:])/50000+e,'b-')\n",
    "ax.plot(xymax,ymax, color='black',linestyle=\"--\",marker=\"o\", markersize=20, linewidth=5)\n",
    "plt.ylabel('Max intensity', fontsize=60)\n",
    "plt.xlabel('Pearls to edge', fontsize=60)\n",
    "plt.xlim(-0.5,27)\n",
    "#plt.ylim(0,11)\n",
    "#plt.grid()\n",
    "plt.xticks(fontsize=45)\n",
    "plt.yticks(fontsize=45)\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(0.005))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "33f128e3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "250"
      ]
     },
     "execution_count": 133,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.convolve(x, np.ones(N)/N, mode='valid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "id": "e709b08d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "125000"
      ]
     },
     "execution_count": 142,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(np.convolve(yvalue, np.ones(10)/10, mode='same'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "id": "7a27f000",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "125000"
      ]
     },
     "execution_count": 136,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(yvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "120c042f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.17322106192806946"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ixi,ixf=find_nearest(np.real(X), 0.28),find_nearest(np.real(X), 0.295)\n",
    "np.max(yvalue[ixi:ixf])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ae71716",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4198cec5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# PLOT of FFT[e] of different pairs of pearls\n",
    "fig,ax = plt.subplots(figsize=(8,10)) # Figure size\n",
    "dataplot=psdfreeM\n",
    "X=dataplot[0]/(33.35641)#*np.sqrt(0.4184))\n",
    "Y=dataplot[1:]\n",
    "lenY=len(Y)-19\n",
    "ixi,ixf=find_nearest(np.real(X), 0.28),find_nearest(np.real(X), 0.295)\n",
    "ymax=[]\n",
    "xymax=[]\n",
    "#plt.figure(figsize=(10,15))\n",
    "#plt.title(\"PSD away from boundary\")\n",
    "colors = iter(cm.jet(np.linspace(0, 1, lenY)))\n",
    "for e in range(0,lenY,1):\n",
    "    yvalue=0.01*np.real(Y[e])/np.average(np.real(Y))\n",
    "    ymax.append(np.max(yvalue[ixi:ixf])*0.5+7.5)\n",
    "    xymax.append(0.35+0.0027*e)\n",
    "    #col = next(colors)\n",
    "    col = ['black','blue','tab:cyan','green','yellow','orange','red']\n",
    "    ax.plot(np.real(X), yvalue-e+6,c=col[e])#, label=\"Ang = {:.1f}\".format(angles[a]))#,colors[e])\n",
    "#for e in range(len(Y)-15):\n",
    "#    plt.plot(np.real(X), 0.15*np.real(Y[e])/np.average(np.real(Y))+e,'b-')\n",
    "    #plt.plot(np.real(ffk[0][len(ffk[0])//2:]), np.real(ffk[e][len(ffk[0])//2:])/50000+e,'b-')\n",
    "#ax.plot(xymax,ymax, color='tab:red')#,linestyle=\"--\")\n",
    "plt.ylabel('PSD (Intensity)', fontsize=20)\n",
    "plt.xlabel('Frequency (cm$^{-1}$)', fontsize=20)\n",
    "plt.xlim(0.15,0.40)\n",
    "plt.ylim(0,11)\n",
    "plt.grid()\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "#ax.xaxis.set_major_locator(ticker.MultipleLocator(0.005))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "167cc751",
   "metadata": {},
   "outputs": [],
   "source": [
    "lll=[]\n",
    "lll.append(np.max(yvalue[ixi:ixf]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "id": "4a70148e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2]"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lll"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d4a7444",
   "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": 5
}
