ELViM_Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89c3bf0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "##Plotting the ELViM projection##\n",
    "\n",
    "##Simple plot with a single color and the position of each structure##\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt ##Imports the libraries needed##\n",
    "from scipy import stats\n",
    "\n",
    "proj1 = 'proj.out' ##Stores the name of projection file##\n",
    "[x1, y1]=np.genfromtxt(proj1, unpack='True') ##Stores the x and y of the projection##\n",
    "\n",
    "fig, ax = plt.subplots() ##Creates the figure##\n",
    "plt.gca().set_aspect('equal') ##Set the aspect of the axis scaling##\n",
    "\n",
    "plot1=ax.scatter(x1, y1, s=3.0, marker='o', c='gainsboro', edgecolor='k', linewidth=0)\n",
    "##Create a scatter plot. One may change the desired color with c='desired_color'##\n",
    "\n",
    "plt.axis('off') ##Removes the x and y axis##\n",
    "plt.savefig('ELViM_projection.png', format='png', dpi=600, bbox_inches = 'tight', pad_inches = 0.05)\n",
    "##Saves the figure file##\n",
    "\n",
    "plt.show() ##Plot the figure##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66a7fe14",
   "metadata": {},
   "outputs": [],
   "source": [
    "##Plotting the ELViM projection with an corresponding information##\n",
    " \n",
    "##ELViM projection with a coordinate associated with each point##\n",
    "##For this script one must add a third column containing the coordinate in the projection file##\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt ##Imports the libraries needed##\n",
    "from scipy import stats\n",
    "\n",
    "projXYZ = 'projection_coordinate.dat' ##Stores the name of the projection file with coordinate## \n",
    "[X, Y, Z]=np.genfromtxt(projXYZ, unpack='True') ##Stores the projection and the coordinate##\n",
    "\n",
    "plotXYZ=ax.scatter(X, Y, s=3.0, marker='o', c=Z, edgecolor='k', linewidth=0, cmap='jet')\n",
    "##Create a scatter plot in which the color corresponds to the coordinate. One may change the name of cmap to\n",
    "#change the colors##\n",
    "\n",
    "plt.axis('off') ##Removes the x and y axis##\n",
    "cbar=fig.colorbar(plotXYZ,ax=ax, pad=0.05, shrink=1) ##Creates a colorbar##\n",
    "cbar.set_label('Density', fontsize=16) ##Set the label of colorbar##\n",
    "cbar.ax.tick_params(labelsize=16) ##Set the size of label##\n",
    "plt.savefig('Density_ab40_monomers.png', format='png', dpi=600, bbox_inches = 'tight', pad_inches = 0.05)\n",
    "##Saves the figure file##\n",
    "\n",
    "plt.show() ##Plot the figure##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cd4ec5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "##Kernel Density Estimation (KDE)##\n",
    "\n",
    "##The KDE can also be calculated for a specific set of points instead of the whole projection##\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt ##Import the libraries##\n",
    "from scipy import stats\n",
    "\n",
    "proj = 'projection.out' ##Save the name of the projection file##\n",
    "x, y = np.genfromtxt(proj, unpack='True') ##Save the coordinates X and Y from the projection file##\n",
    "xy=np.vstack([x,y]) ##Concatenate the sequences vertically##\n",
    "\n",
    "z=stats.gaussian_kde(xy, weights=,bw_method=0.03)(xy) ##Estimate the probability density function of the points## \n",
    "##One should change the bw_method accordingly with the number os points in the projection##\n",
    "\n",
    "idx=np.argsort(z) ##Creates an array of indices##\n",
    "\n",
    "fig, ax = plt.subplots() ##Create the figure and a set of subplots##\n",
    "plt.gca().set_aspect('equal') ##Set the aspect of the axis scaling##\n",
    "sc=ax.scatter(x[idx], y[idx], s=3.0, marker='o', c=np.log(z[idx]),edgecolor='k', linewidth=0.02, alpha=1 , cmap='afmhot_r')\n",
    "##Creates a scatter plot in which the color corresponds to local density##\n",
    "\n",
    "plt.axis('off') ##Removes the x and y axis##\n",
    "cbar=fig.colorbar(sc,ax=ax, pad=0.05, shrink=1) ##Creates a colorbar##\n",
    "cbar.set_label('Density', fontsize=16) ##Set the name of the colorbar##\n",
    "cbar.ax.tick_params(labelsize=16) ##Set the size of the name of the colorbar##\n",
    "plt.savefig('kde_ln.png', format='png', dpi=600, bbox_inches = 'tight', pad_inches = 0.05) ##Save the figure##\n",
    "plt.show() ##Plot the figure##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3226a3d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "##Entropy Estimation##\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt ##Import the libraries##\n",
    "from scipy import stats\n",
    "\n",
    "proj = 'projection.out' ##Save the name of the projection file##\n",
    "x, y, = np.genfromtxt(proj, unpack='True') ##Save the coordinates X and Y from the projection file##\n",
    "\n",
    "w=np.ones(len(x)) ##Assigns the same weight to all points##\n",
    "stat,x_edge, y_edge,  binnumber = stats.binned_statistic_2d(x, y, w, 'count', bins=[25,25],expand_binnumbers=True)\n",
    "##\n",
    "\n",
    "############  countagem de amostras por bin para tipo 1\n",
    "w=np.ones(len(x)) ##Assigns the same weight to all points##\n",
    "stat1,_, _ , binn1 = stats.binned_statistic_2d(x, y, w, 'count', bins=[x_edge, y_edge],expand_binnumbers=True)\n",
    "##Set the grid and count the number of points in each bin##\n",
    "\n",
    "stat1=stat1/np.sum(stat1) ##Normalizes the function##\n",
    "stat1[stat1 == 0] = 1 ##Assigns all zeros equal 1 to compute de log##\n",
    "H = -np.sum(stat1*np.log(stat1)) ##Compute de log of the distribuction##\n",
    "print(H) ##Return the value of the entropy on the screen##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8ada089",
   "metadata": {},
   "outputs": [],
   "source": [
    "##Comparison between differents KDE##\n",
    "\n",
    "##The KDE can be calculated for a specific set of points instead of the whole projection##\n",
    "##With this analyse one can compare two differents sets of KDE that came from the same projection##\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import seaborn as sns; sns.set(style=\"white\",color_codes=True)\n",
    "\n",
    "xall, yall = np.genfromtxt('proj_all.dat', unpack=True) ##Loads the whole projection to define single grid##\n",
    "x1, y1 = np.genfromtxt('ab40_tetramers.dat', unpack=True) ##Loads the set of the first KDE##\n",
    "x2, y2 = np.genfromtxt('ab42_tetramers.dat', unpack=True) ##Loads the set of the second KDE##\n",
    "values=np.vstack([xall,yall]) ##Concatenate the whole projection##\n",
    "\n",
    "xmax=xall.max()+0.05            \n",
    "xmin=xall.min()-0.05     ##Defines the grid maxima and minimus##\n",
    "ymax=yall.max()+0.05     ##The 0.05 ensures that there is no edge effect\n",
    "ymin=yall.min()-0.05\n",
    "\n",
    "X,Y = np.mgrid[xmin:xmax:300j, ymin:ymax:300j]  ##Set the grid with (300 x 300) bins\n",
    "positions = np.vstack([X.ravel(), Y.ravel()])   ##Concatenates the X and Y position in a scipy format##\n",
    "\n",
    "values1 = np.vstack([x1,y1])                    \n",
    "kernel1 = stats.gaussian_kde(values1) ##Computes the KDE for the first set##\n",
    "kernel1.set_bandwidth(bw_method=kernel1.factor*0.4)\n",
    "c1=kernel1.evaluate(values)                        ##KDE value of each point##\n",
    "Z1 = np.reshape(kernel1(positions).T, X.shape)     ##KDE value in each bin##\n",
    "\n",
    "values2 = np.vstack([x2,y2])\n",
    "kernel2 = stats.gaussian_kde(values2)\n",
    "kernel2.set_bandwidth(bw_method=kernel2.factor*0.4) ##Computes the KDE for the second set##\n",
    "c2=kernel2.evaluate(values)\n",
    "Z2 = np.reshape(kernel2(positions).T, X.shape)\n",
    "\n",
    "Z=Z1-Z2 ##Computes the difference between the KDE##\n",
    "lim=max(abs(np.max(Z)), abs(np.min(Z))) ##Applies a limit to make the scale symmetrical##\n",
    "plt.figure() ##Creates the figure##\n",
    "plt.scatter(xall, yall, s=1.0, marker='o',c='k', linewidth=0, alpha=0.05) ##Plot the points of both sets for reference##\n",
    "img=plt.imshow(Z.T, origin='lower',extent=([xmin, xmax, ymin, ymax]), cmap=\"RdBu\") ##Plot the KDE differences##\n",
    "plt.clim(-lim, lim)   ##Defines the zero as the center of the scale##\n",
    "plt.colorbar(img)     ##Creates a colorbar##\n",
    "plt.axis('off')       ##Remove the X and Y axis##\n",
    "plt.show()            ##Plot the whole figure##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4c27d3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "##Free energy estimation##\n",
    "\n",
    "##To compute the free energy estimation from the elvim projection some thermodynamic analysis are necessary\n",
    "#We recommend the Weighted Histogram Analysis Method (WHAM) to obtain the necessary file\n",
    "#A tutorial on how to execute WHAM can be found in the SMOG 2 manual.##\n",
    "\n",
    "\n",
    "##After running WHAM, you will get a \"dos\" file, which contains the potential energy, a reaction coordinate and\n",
    "#the density of states.\n",
    "\n",
    "##To run this script one must provide the projection file, a file containing the potential energy of each structure\n",
    "#in the projection and the \"dos\" file from WHAM##\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "\n",
    "kb=0.0019872041 ##Boltzmann constant##\n",
    "se = np.loadtxt('energy.dat') ##Stores the potential energy of each structure##\n",
    "e, q, dos =np.genfromtxt('dos', unpack='True') ##Stores the energy, coordinate and the density of states##\n",
    "\n",
    "ev=np.unique(e) ##Set the left boundary of energy bins##\n",
    "qv=np.unique(q) ##Set the left boundary of coordinate bins##\n",
    "neb=len(ev)     ##Set the number of energy bins##                \n",
    "nqb=len(qv)     ##Set the number of coordinate bins## \n",
    "ddos=np.zeros([neb, nqb]) #Set a matrix of zeros to store the bins##\n",
    "\n",
    "for n in range(len(dos)):        #Rewrite the dos file in a matrix format##\n",
    "    i, j =np.where(ev == e[n] )[0], np.where(qv == q[n])[0]\n",
    "    ddos[i,j] = dos[n]\n",
    "\n",
    "dosE = np.sum(ddos, axis=1)   ##Sum the coordinate bins, leaving only the energy bins##\n",
    "dosT=np.sum(dosE)             ## dos total ##\n",
    "\n",
    "ie = np.digitize (se, ev)  ##Set the number of energy bins in each frame (starting at 1)\n",
    "s = len(se)                ##Set the number of points in the projection##     \n",
    "w=np.zeros(s)              ##Declares a variable for the weight of each point##\n",
    "    \n",
    "for i in range(s):\n",
    "    w[i]=dosE[ie[i]-1]/dosT*np.exp(-se[i]/(kb*300)) ##Compute de weight of each point##\n",
    "w=w/(np.sum(w))                                    ##Normalize the value##\n",
    "\n",
    "x, y =np.genfromtxt('projection_linear2.out', unpack=True) ##Loads the projection file##\n",
    "xy=np.vstack([x,y]) ##Concatenates the values##\n",
    "z=stats.gaussian_kde(xy, weights=w)(xy) ##Computes the density estimation##\n",
    "f=-kb*300*np.log(z) ##Computes the free energy##\n",
    "idx=np.argsort(f)\n",
    "\n",
    "fig, ax = plt.subplots() ##Creates the figure##\n",
    "plt.gca().set_aspect('equal') ##Set the aspect of the axis scaling##\n",
    "sc=ax.scatter(x[idx], y[idx], s=3.0, marker='o', c=f[idx], edgecolor='k', linewidth=0.02, alpha=1 , cmap='Blues')\n",
    "##Creates a scatter plot in which the color corresponds to free energy##\n",
    "\n",
    "plt.axis('off') ##Removes the x and y axis##\n",
    "cbar=fig.colorbar(sc,ax=ax, pad=0.05, shrink=1) ##Set a colorbar##\n",
    "cbar.set_label('Free_projection', fontsize=16) ##Set the name of the colorbar##\n",
    "cbar.ax.tick_params(labelsize=16) ##Set the size of the label##\n",
    "plt.savefig('kde_ln.png', format='png', dpi=600, bbox_inches = 'tight', pad_inches = 0.05) ##Saves the figure##\n",
    "plt.show() ##Plot the figure##"
   ]
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}