{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "160a91f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.fft import fft, ifft, fftn, ifftn\n",
    "import scipy\n",
    "import glob\n",
    "import sys\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.collections import LineCollection\n",
    "from matplotlib import cm\n",
    "#from IPython.display import Image\n",
    "#from mpl_toolkits.mplot3d import Axes3D\n",
    "#from functools import partial\n",
    "#import pathos.multiprocessing as mp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "24eb0df0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# are the interaction energies of different pairs dabco-FeP (see trajectories for details)\n",
    "sep1=np.loadtxt(\"../gromacs/pdp2x14x/separa/minienergy3.pdp.28x.yy.dat\")\n",
    "sep2=np.loadtxt(\"../gromacs/pdp2x14x/separa/minienergy4.pdp.28x.yy.dat\")\n",
    "sep3=np.loadtxt(\"../gromacs/pdp2x14x/separa/minienergy5.pdp.28x.yy.dat\")\n",
    "sep4=np.loadtxt(\"../gromacs/pdp2x14x/separa/minienergy6.pdp.28x.yy.dat\")\n",
    "sep=np.vstack((sep1,sep2,sep3,sep4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "4fde6809",
   "metadata": {},
   "outputs": [],
   "source": [
    "# minienergy7 and minienergy8 are the interaction energies between a dabco molecule and the FeP site\n",
    "sep1=np.loadtxt(\"../gromacs/pdp2x14x/separa/minienergy7.pdp.28x.yy.dat\")\n",
    "sep2=np.loadtxt(\"../gromacs/pdp2x14x/separa/minienergy8.pdp.28x.yy.dat\")\n",
    "sep=np.vstack((sep1,sep2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "ec077a2f",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Separation ($\\\\AA$)')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 61712 characters]\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Overview surface and atomic indexes\n",
    "fig,ax = plt.subplots(figsize=(8,8)) # Figure size\n",
    "X=np.asarray(list(range(0, 50)))/10.\n",
    "range(0,50)\n",
    "#plt.title('Eigenmodes for {:d} Pearls: -{:d}-[-{:d}-]-{:d}- '.format(Np,k2,k1,k2))\n",
    "#for e in range(0,5):\n",
    "#plt.plot(X,np.real(Y[e]),'k--')\n",
    "colors=['g','r','m','b','k','y','r','c','g','m'] #mbcgyrk\n",
    "for e in range(len(sep)):\n",
    "    Y=sep[e]-sep[0][0]\n",
    "    ax.plot(X,Y,c=colors[e],label=\"Separation \"+str(e+1))\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "#plt.ylim(70,100)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=16)\n",
    "plt.ylabel('Potential Energy (kJ/mol)', fontsize=20)\n",
    "plt.xlabel('Separation ($\\AA$)', fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f97aece1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "29.092894999999864\n",
      "23.055299999999896\n",
      "31.757394999999907\n",
      "57.1036349999999\n"
     ]
    }
   ],
   "source": [
    "ene=[]\n",
    "for e in range(len(sep)):\n",
    "    Y=sep[e]-sep[0][0]\n",
    "    aveY=np.average(Y[30:])\n",
    "    ene.append(aveY)\n",
    "    print(aveY)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9273d051",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(14.726764999999816, 10.993035000000106, 8.32853500000007, 6.037594999999968)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# it seems something is wrong here\n",
    "A,B,C,D=ene[0],ene[1],ene[2],ene[3]\n",
    "F2=D-2*B\n",
    "F4=A-B\n",
    "F1=C+3*B-A-D\n",
    "F3=A+D-2*B-C\n",
    "F1,F2,F3,F4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1425e638",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10.993035000000106, 14.726764999999816, 6.037594999999968, 8.32853500000007)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# done on 30.03.23\n",
    "# A = F2+F3+F4\n",
    "# B = F2+F4\n",
    "# C = F1+F2+F3\n",
    "# D = F1+2*F2+2*F4\n",
    "######\n",
    "A,B,C,D=ene[0],ene[1],ene[2],ene[3]\n",
    "F1=D-2*B\n",
    "F3=A-B\n",
    "F2=C+3*B-A-D\n",
    "F4=A+D-2*B-C\n",
    "F1,F2,F3,F4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c208a0ff",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2945.3529999999632, 2198.607000000021, 1665.707000000014, 1207.5189999999934)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F1*200,F2*200,F3*200,F4*200"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "id": "b66ee9d4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(39.02683767294731, 29.132222350124948, 22.071132628141143)"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F1/F4*16, F2/F4*16, F3/F4*16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfedbaf2",
   "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
}
