tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e2d1b9bd-a8a0-4cf1-95a6-7368eebbb0ec",
   "metadata": {},
   "source": [
    "# What is metadynminer"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e34e2aae-ab72-491a-9113-dadd3ddac974",
   "metadata": {},
   "source": [
    "Welcome to Metadynminer, our Python package developed to make analysis of metadynamics simulations easy and user-friendly. This Jupyter notebook is designed to introduce metadynminer to new users and also to explain the algorithms used by metadynminer to some extent. \n",
    "\n",
    "Just a reminder, if you are using this notebook in our online service to process your own results, do not forget to download all files you may have created to your PC for later use. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d591434-dccf-4632-b421-b04fb13ba0c8",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import metadynminer as mm\n",
    "print(f\"Loaded metadynminer version {mm.__version__}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ab202cd-a0f8-4c5b-b66e-bdbdf8ab2e3d",
   "metadata": {},
   "source": [
    "You can uncomment and run the next line to enable ```%matplotlib widget``` for better interactivity with the graphs, if you want to. It is not necessary in most cases. To turn it back off, you may need to restart the kernel of this notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bff99bc-f1ca-4e80-b620-50a47b4aadb8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %matplotlib widget"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04c50bc0-7644-4355-941f-4fba7ec2f501",
   "metadata": {},
   "source": [
    "# Load your HILLS file "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73cb4b76-7bcc-4bfd-9dc4-b4a67a36bd4b",
   "metadata": {},
   "source": [
    "First upload your HILLS file and load it with the line below. If the name of your HILLS file is not the default \"HILLS\", specify it to the ```name``` keyword. \n",
    "\n",
    "If you do not have your own HILLS file, you can try our code with our HILLS file from the oxytocin metadynamics simulation used for testing and demonstration, using \n",
    "```python\n",
    "hills = mm.Hills(name=\"oxytocin\")\n",
    "```\n",
    "\n",
    "Metadynminer automatically detects which CVs are periodic, however, if it fails to detect periodicity correctly, \n",
    "you can specify the periodicity of your CVs as boolean Python list\n",
    "(for example, if you have two periodic CVs, set ```periodic = [True, True]```): \n",
    "```python\n",
    "hills = mm.Hills(name=\"HILLS\", periodic=[True, True])\n",
    "```\n",
    "\n",
    "Default periodicity of periodic CVs is set from -$\\pi$ to $\\pi$. If your CVs have different periodicity, you have to specify it to the respective keyword ```cv1per```, ```cv2per``` or ```cv3per```. \n",
    "\n",
    "An example for CV 1 with periodicity from 0 to 2$\\pi$: \n",
    "```python\n",
    "import numpy as np\n",
    "hills = mm.Hills(name=\"HILLS\", periodic=[True, True], cv1per=[0, 2*np.pi])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62a89391-18e5-4203-86fd-3eda3a238d80",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "hills = mm.Hills(name=\"oxytocin\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c89d939f-ecfd-45aa-9142-e50fd2a374b2",
   "metadata": {},
   "source": [
    "# Plot CVs and heights of the hills"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39189356-3894-47b6-ac80-ea666e886c2a",
   "metadata": {},
   "source": [
    "You can show the evolution of your molecular system by plotting values of the collective variables in time. For this, use the ```plot_CV``` method. Specify which CV should be plotted to the ```CV``` keyword. If you want to save and download the image with the plot later, specify the name for the file to the  ```png_name``` keyword, for example: \n",
    "```python\n",
    "hills.plot_CV(png_name=\"CV_1_plot.png\", CV=1)\n",
    "```\n",
    "Axis labels will be generated automatically based on the information inside the HILLS file. If you want to use your own labels, set them as strings to ```xlabel```, ```ylabel``` keywords:\n",
    "```python\n",
    "hills.plot_CV(png_name=\"CV_1_plot.png\", CV=1, ylabel=\"CV 1\")\n",
    "```\n",
    "You can also change the time unit to be shown as necessary: ```tu=\"ns\"``` etc. \n",
    "\n",
    "To export the image in publication quality, you might want to use other keywords to set the size and dpi of the image, for example: ```dpi=300```, ```image_size=[16,9]``` and ```image_size_unit=\"cm\"```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58969d64-84ef-4888-9481-454cdd000e89",
   "metadata": {},
   "outputs": [],
   "source": [
    "hills.plot_CV(png_name=None, CV=1, tu=\"ns\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "00d701a3-c188-4eec-9a53-6a82b52454ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "hills.plot_CV(png_name=None, CV=2, tu=\"ns\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba2bb070-c029-4a94-b383-d55e3f06d02b",
   "metadata": {},
   "source": [
    "In well-tempered metadynamics, it might be also usefull to visualise the heights of the hills added during the simulation, \n",
    "as these tend to decrease over time when the system revisits already explored parts of the CV space: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82672dcf-c251-4640-bc2f-d917888cd6a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "hills.plot_heights(png_name=None, energy_unit=\"kJ/mol\", title=\"heights of hills\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13f775f1-935a-4bb0-9748-c012b64d566b",
   "metadata": {},
   "source": [
    "# Calculate the free energy surface from hills file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a44a4c35-9205-46b3-8c94-70356daf6d10",
   "metadata": {},
   "source": [
    "You have two options to choose from here: \n",
    "If you use ```original=False```, the FES will be calculted by the Hillsum algorithm - only one Gaussian hill will be precalculated \n",
    "and then moved and scaled and added to the FES. This is less precise, but much faster than calculating each Gaussian explicitly. \n",
    "This is ideal for casual visualisation purposes. \n",
    "\n",
    "If you use ```original=True```, each Gaussian hill will be calculated explicitly, which is more computationally demanding, but also exact. \n",
    "This method is giving the same results as the ```sum_hills``` function in Plumed (tested with Plumed v2.8.0)\n",
    "\n",
    "If you don't need high resolution FES and/or you find the calculation too slow, you can decrease the ```resolution``` of the FES. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7006d903-9804-4bc9-8c8d-819cc96f4037",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = mm.Fes(hills, original=False, resolution=256)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd05470c-cf8c-4b3a-8e90-a3a8965eed19",
   "metadata": {},
   "source": [
    "Then you can visualise the FES. Keywords controlling the image size, dpi, name etc. are the same as for the CV plots. \n",
    "\n",
    "For 2D and 3D FESs, you can change the colormap to any of the matplotlib colormaps, for example ```cmap=\"rainbow\"```. \n",
    "\n",
    "On 2D and 3D FESs, you can also control the contours or isosurfaces to be shown, either as a list specifying the values of free energy to be visualised as here: ```levels=[20,40,60]```, or by specifying the spacing between each contour/isosurface: ```contours_spacing=20```. \n",
    "\n",
    "Sometimes it is useful to change the minimum or maximum value of free energy to be colored with keywords ```vmin```, ```vmax```.\n",
    "\n",
    "To plot only a slice of the FES, you can use keywords ```xlim=[x_minimum, x_maximum]```, ```ylim=[y_minimum, y_maximum]``` with a list of the beginning and the end of the slice.\n",
    "\n",
    "For 3D FESs, it may be useful to change the opacity of the isosurfaces to get optimal visibility: ```opacity=0.2```. \n",
    "\n",
    "\n",
    "For example:\n",
    "```python\n",
    "fes.plot(png_name=None, contours=True, cmap = \"RdYlBu_r\", \n",
    "             energy_unit=\"kJ/mol\", xlabel=None, ylabel=None, zlabel=None, label_size=12, clabel_size = 12,\n",
    "             image_size=[9,6], image_size_unit=\"in\", dpi=100, vmin = 0, vmax = None, levels=[20,40,60,80,100])\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33a11a5b-8a8f-48fa-b29e-220c6758495e",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes.plot( cmap = \"RdYlBu_r\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d9e2b01-433c-4903-ace7-5610d3875541",
   "metadata": {},
   "source": [
    "# Identify free energy minima"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5a10071-5a58-4cfc-b29d-79d423b60341",
   "metadata": {},
   "source": [
    "You have two options: \n",
    "If you set the keyword ```precise=True```, the local minima will be identified by an algorithm which: \n",
    "1. finds all local minima, even very shallow and probably unimportant minima, \n",
    "2. each point on the FES will be assigned to the minimum the system would most likely go to, \n",
    "   if it only follows the gradient of free energy, and \n",
    "3. free energy value of minima will be calculated from each point on FES assigned to the respective minima. \n",
    "   This results in more precise free energy values, as it accounts for the width of the minimum as well. \n",
    "   For this calculation the unit of free energy (```energy_unit=\"kJ/mol\"``` or ```energy_unit=\"kcal/mol\"```) and \n",
    "   the thermodynamical temperature (```temp```) of the simulation must be supplied. \n",
    "   This algorithm does not use the ```nbins``` keyword. \n",
    "\n",
    "Example:\n",
    "```python\n",
    "minima = mm.Minima(fes, precise=True, temp=300.0, energy_unit=\"kJ/mol\")\n",
    "```\n",
    "\n",
    "If you set ```precise = False```, the method will use the original algorithm from the metadynminer package for R. \n",
    "In this algorithm the FES is first divided to number of bins (can be set with option ```nbins```, default is 8 for 1D FES, $8\\times 8$ for 2D FES or $8\\times 8 \\times 8$ for 3D FES)\n",
    "and the absolute minima is found for each bin. Then the algorithm checks \n",
    "if this point is really a local minimum by comparing to the surrounding points of FES.\n",
    "This algorithm only accounts for the depth of each minima, which is less precise, but usually sufficient. \n",
    "\n",
    "In some cases, this algorithm is the prefered one, \n",
    "because on some free energy landscapes the total number of local free energy minima can reach tens of thousands, \n",
    "which makes the calculation using precise algorithm slow and impractical. \n",
    "\n",
    "Example: \n",
    "```python\n",
    "minima = mm.Minima(fes, precise=False, nbins=8)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6139d9e3-5b94-4417-ad72-22f7b96def03",
   "metadata": {},
   "outputs": [],
   "source": [
    "minima = mm.Minima(fes, precise=True, temp=300.0, energy_unit=\"kJ/mol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dee99b64-a657-4f04-9317-87ae5260fadb",
   "metadata": {},
   "source": [
    "Print the list of the local minima:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93f83df5-5d4b-4396-bc49-826cc3f58ae4",
   "metadata": {},
   "outputs": [],
   "source": [
    "minima.minima"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fc4c016-8444-47f6-aeaa-384008419f25",
   "metadata": {},
   "source": [
    "You can visualise the FES with the local minima as letters. All keywords are analogous as for the ```fes.plot()``` function.\n",
    "The color of the letters changes automatically to ensure their good visibility, but if you want to override this behaviour, provide a matplotlib color like this: ```color=\"black\"```. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "969b2add-4f83-410d-84e9-fc750e308a1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the free energy surface with minima\n",
    "minima.plot(contours_spacing=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d12bea4-76c6-4474-bcd7-e8f6bd55f336",
   "metadata": {},
   "source": [
    "# Construct free energy profile"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b096e969-45e7-4190-a001-cc822454c4b9",
   "metadata": {},
   "source": [
    "In metadynamics simulations, it is useful to visualise how the FES was changing during the simulation, which can suggest (not prove) whether the FES is converged or not and the simulation should be prolonged. \n",
    "Here, you should provide the list of minima as well as the hills objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "355f6d48-0517-421f-9d87-435873df516c",
   "metadata": {},
   "outputs": [],
   "source": [
    "prof = mm.FEProfile(minima,hills)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4f93cee-67e1-4d29-a98a-76bd2d318d80",
   "metadata": {},
   "source": [
    "Plot the free energy profile:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fed8bfe-fa73-48f0-a04f-2f7ac5c525d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "prof.plot(legend=True, png_name=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a263549f-72c9-4a91-971e-52e49183cd43",
   "metadata": {},
   "source": [
    "# Other optional methods for analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4ac1ef7-050d-49c4-964e-49fa84284331",
   "metadata": {},
   "source": [
    "You can remove one CV from an existing FES (for example to make visualisation of 3D FES easier). For this purpose use the ```Fes.remove_CV()``` method, which will return a new FES with one CV removed. The algorighm converts the FES to probabilities, then sums the probabilities along the given ```CV``` to be removed and then converts these sums back to free energy values. Because of this, you should provide the temperature of the simulation as well as the unit of free energy used in HILLS file. ```temp=300.0``` Kelvin and ```energy_unit=\"kJ/mol\"``` are the default values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73f38372-bd8e-4157-96ca-8be1cce7ef08",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes_CV1 = fes.remove_CV(CV=2, temp=300.0, energy_unit=\"kJ/mol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0831c47-8ae7-4b6f-8201-0b23fd07b097",
   "metadata": {},
   "source": [
    "You can work with the new fes object in the same way as with the original one, e. g. find minima, plot them etc.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dbc1ae85-16f4-45ef-a193-65f53df2467a",
   "metadata": {},
   "outputs": [],
   "source": [
    "minima_1 = mm.Minima(fes_CV1)\n",
    "minima_1.minima"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c214175-1d06-496c-a2df-1eeb3250aa21",
   "metadata": {},
   "outputs": [],
   "source": [
    "minima_1.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd31b3e3-5c65-46a5-8411-1bbfa3644c47",
   "metadata": {},
   "source": [
    "Another alternative way to visualise 2D FES is by creating a surface plot. This method only works for 2D FESs. This works best together with ```%matplotlib widget``` turned on (you can find the line at the beginning of this notebook). If you find the animation too slow, it may be necessary to decrease the resolution of the FES. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd293e02-f560-47d5-a5f4-ef153528adbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes.surface_plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd9d1469-e79d-4ff4-9527-5f867e165c45",
   "metadata": {},
   "source": [
    "Visualising 3D FES can be challenging. One possible way is by plotting different isosurfaces with specific free energy values in 3D space, as in the ```Fes.plot()``` and ```Minima.plot()``` methods. Another way is by creating an animation showing different isosurfaces at diferent times with ```Fes.make_gif()``` method and it's analog ```Minima.make_gif()```. These methods are only available for 3D FESs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68a44536-ce1e-4a01-8255-a8d41f11ade5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes.make_gif(gif_name=\"fes_animation.gif\", energy_unit=\"kJ/mol\",\n",
    "                  opacity=0.2, frames=64)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbbaf80a-a52a-4ff2-8e4e-ecf4985f5ca0",
   "metadata": {},
   "source": [
    "The resulting animation (if available). Maybe you will need to refresh this page of your browser for the animation to be shown:\n",
    "\n",
    "<img src=\"fes_animation.gif\" width=\"750\" align=\"center\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f678d7a8-f132-4135-9e91-ee50a1e859bb",
   "metadata": {},
   "source": [
    "It is also possible to create an animation of the flooding of the free energy surface as the simulation ran. This can be useful to visualise the convergence of the FES and the evolution of the CV values during the simulation at the same time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d613dd8-88c2-4f86-a896-1ebc3b4a12b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "fes.flooding_animation(step=1000, gif_name=\"flooding.gif\", fps=10, enable_loop=True, \n",
    "                       contours_spacing=10, with_minima=True, use_vmax_from_end=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1c0752e-86ca-41c2-9caa-21c7bbf80fbd",
   "metadata": {},
   "source": [
    "Again, maybe you will need to refresh this page of your browser for the resulting animation to be shown:\n",
    "\n",
    "<img src=\"flooding.gif\" width=\"750\" align=\"center\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7dbb7306-a521-447f-bb71-35f872af44e6",
   "metadata": {},
   "source": [
    "# Reweighting"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cca1df8-baa6-46a7-909f-cf33e80feeeb",
   "metadata": {},
   "source": [
    "There is also more precise (and sometimes more complicated) way to obtain the free energy surface, commonly reffered to as *reweighting*. Here, we show Tiwary's reweighting. In principle, FES can be calculated with respect to any CV by calculating the free energy of states with some value from their relative population compared to other states. That is quite simple for unbiased simulations. In metadynamics, some states are discouraged more than others by the bias potential, which must be accounted for during reweighting. In other words, for states that were populated even though they were discouraged by bias potential, the relative population has higher weight compared to not-so-discouraged states. Reweighting can be used to calculate the FES for the same CVs that were biased by metadynamics, but also with respect to different CVs, as long as you can provide the values of the CVs during simulation. You can calculate these values for example with ```plumed driver```. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43c05470-d03f-4bf8-b226-faeb74d887b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c5977bb-9c7d-43cd-b8c9-112cfffec6d5",
   "metadata": {},
   "source": [
    "Specify temperature of simulation, energy unit in HILLS file and resolution of resulting reweighted FES:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8eee55f9-0af4-4d9a-833e-cf3a9d25348d",
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature = 300.0\n",
    "energy_unit=\"kJ/mol\"\n",
    "resolution=50 \n",
    "maxfes = 75"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6997739c-5a34-4284-81e0-40caee0609f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "if energy_unit == \"kJ/mol\":\n",
    "    kT = 8.314*temperature/1000\n",
    "elif energy_unit == \"kcal/mol\":\n",
    "    kT = 8.314*temperature/1000/4.184\n",
    "biasfactor = float(hills.biasf[0])\n",
    "nsteps=50\n",
    "outfes = mm.Fes(hills, resolution=resolution, calculate_new_fes=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c94a71f7-40e0-48c5-af8c-1e281a97638a",
   "metadata": {},
   "source": [
    "Here, you should provide one other output file from metadynamics simulation (common name being ```COLVAR```) containing the total bias potential present applied to the system at each sampled  moment. If you also calculated the reweighting factor on the fly (with ```CALC_RCT``` keyword), the following calculation will be a bit more simple. The default name of the column containing the reweighting values is ```metad.rbias```. \n",
    "\n",
    "If your ```COLVAR``` file does not contain the ```rbias``` values, but only ```metad.bias```, you should calculate the Tiwary's correction manually - skip to the subsection 8.2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c49261c-d617-44d2-80c7-6bcb58b5c1ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "colvar = np.loadtxt(\"COLVAR_oxytocin\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf95b281-10d0-4679-9c51-3627244c98a4",
   "metadata": {},
   "source": [
    "## Algorithm using the reweighting factor ```metad.rbias``` calculated on-the-fly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b847c389-5468-4bd2-8085-05a2ac03b013",
   "metadata": {},
   "outputs": [],
   "source": [
    "# excluding values where time == 0.0, because otherwise when working with concatenated simulations, the same state would be sampled twice \n",
    "# (once at the end of previous simulation part and then at the begining of the continuation)\n",
    "colvar1 = colvar[colvar[:,0]!=0.0,1] # colvar[:,1] contains values of CV1 during simulation\n",
    "colvar2 = colvar[colvar[:,0]!=0.0,2] # colvar[:,2] contains values of CV2 during simulation\n",
    "rbias = colvar[colvar[:,0]!=0.0,4] # colvar[:,3] contains the metad.rbias column\n",
    "\n",
    "nsamples = colvar[colvar[:,0]!=0.0].shape[0]\n",
    "step = np.arange(0, nsamples)*nsteps/(nsamples+1)\n",
    "ix = (resolution*(colvar1-outfes.cv1min)/(outfes.cv1max-outfes.cv1min)).astype(int)\n",
    "iy = (resolution*(colvar2-outfes.cv2min)/(outfes.cv2max-outfes.cv2min)).astype(int)\n",
    "probabilities = np.zeros((outfes.fes.shape))\n",
    "ebias = np.exp(rbias/kT)\n",
    "for i in range(int(nsamples)):\n",
    "    probabilities[ix[i], iy[i]] = probabilities[ix[i],iy[i]] + ebias[i]\n",
    "\n",
    "probabilities[probabilities == 0.] = np.nextafter(0.,1)\n",
    "outfes.fes = -kT*np.log(probabilities)\n",
    "outfes.fes = outfes.fes - np.min(outfes.fes)\n",
    "outfes.fes[outfes.fes>maxfes] = maxfes\n",
    "outfes.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21702bf2-d4e3-4414-99fe-f001e2d815e9",
   "metadata": {},
   "source": [
    "## Algorithm using ```metad.bias``` calculating Tiwary's correction post *ex post*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8318a66f-769f-4850-b6e3-8d2d268cd3c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "hills.hills.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1d1934f-7bf2-49bc-9e7a-633915fdf360",
   "metadata": {},
   "outputs": [],
   "source": [
    "step = np.arange(1,nsteps+1)*hills.cv1.shape[0]/nsteps\n",
    "s1 = np.zeros((step.shape))\n",
    "for i in range(len(s1)):\n",
    "    # print progress \n",
    "    print(f\"Progress: {((i+1)/len(s1)):.2%} finished. \", end=\"\\r\")\n",
    "    if i == 0:\n",
    "        s1[i] = np.sum(np.exp((-mm.Fes(hills, original=False, print_output=False, \n",
    "                                   time_max=step[i], subtract_min=False).fes/kT)*((biasfactor-1)/biasfactor)))\n",
    "    else:\n",
    "        s1[i] = s1[i-1] + np.sum(np.exp((-mm.Fes(hills, original=False, print_output=False, \n",
    "                                   time_min=step[i-1], time_max=step[i], subtract_min=False).fes/kT)*((biasfactor-1)/biasfactor)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdbcda49-da23-4685-8c2c-7ef42fc17905",
   "metadata": {},
   "outputs": [],
   "source": [
    "s2 = np.zeros((step.shape))\n",
    "for i in range(len(s2)):\n",
    "    # print progress \n",
    "    print(f\"Progress: {((i+1)/len(s1)):.2%} finished. \", end=\"\\r\")\n",
    "    if i == 0:\n",
    "        s2[i] = np.sum(np.exp((-mm.Fes(hills, original=False, print_output=False, \n",
    "                                   time_max=step[i], subtract_min=False).fes/kT/biasfactor)*((biasfactor-1)/biasfactor)))\n",
    "    else:\n",
    "        s2[i] = s2[i-1] + np.sum(np.exp((-mm.Fes(hills, original=False, print_output=False, \n",
    "                                   time_min=step[i-1], time_max=step[i], subtract_min=False).fes/kT/biasfactor)*((biasfactor-1)/biasfactor)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "115bd001-0d6c-48f4-82dc-7b45c1fc890f",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ebetac = s1/s2\n",
    "\n",
    "# excluding values where time == 0.0, because otherwise when working with concatenated simulations, the same state would be sampled twice \n",
    "# (once at the end of previous simulation part and then at the begining of the continuation)\n",
    "colvar1 = colvar[colvar[:,0]!=0.0,1] # colvar[:,1] contains values of CV1 during simulation\n",
    "colvar2 = colvar[colvar[:,0]!=0.0,2] # colvar[:,2] contains values of CV2 during simulation\n",
    "bias = colvar[colvar[:,0]!=0.0,3] # colvar[:,3] contains the metad.bias column\n",
    "\n",
    "nsamples = colvar[colvar[:,0]!=0.0].shape[0]\n",
    "step = (np.arange(0, nsamples)*nsteps/(nsamples+1)).astype(int)\n",
    "ix = (resolution*(colvar1-outfes.cv1min)/(outfes.cv1max-outfes.cv1min)).astype(int)\n",
    "iy = (resolution*(colvar2-outfes.cv2min)/(outfes.cv2max-outfes.cv2min)).astype(int)\n",
    "ebetac = np.repeat(ebetac, int(bias.shape[0]/ebetac.shape[0]))\n",
    "ebias = np.exp(bias/kT)/ebetac\n",
    "probabilities = np.zeros((outfes.fes.shape))\n",
    "for i in range(int(nsamples)):\n",
    "    probabilities[ix[i], iy[i]] = probabilities[ix[i],iy[i]] + ebias[i]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3423ce9-c28b-46ab-a8f0-6953a657248f",
   "metadata": {},
   "outputs": [],
   "source": [
    "probabilities[probabilities == 0.] = np.nextafter(0.,1)\n",
    "outfes.fes = -kT*np.log(probabilities)\n",
    "outfes.fes = outfes.fes - np.min(outfes.fes)\n",
    "outfes.fes[outfes.fes>maxfes] = maxfes\n",
    "outfes.plot()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "base"
  },
  "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.9"
  },
  "toc-autonumbering": true,
  "toc-showcode": false,
  "toc-showmarkdowntxt": true
 },
 "nbformat": 4,
 "nbformat_minor": 5
}