05_analysis/01_internal_coordinates_averages.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Internal coordinates and averages\n",
    "\n",
    "This notebooks focuses on internal coordinates, distances, angle and dihedral angles. The name \"internal coordinates\" refers to the fact that they are insensitive to global rotations or translations of  the system. \n",
    "\n",
    "These internal coordinates will also be used to demonstrate the computation of averages of time series and the error on such averages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "before-import"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.patches as mpatches\n",
    "import matplotlib.pyplot as plt\n",
    "import mdtraj\n",
    "import nglview\n",
    "import numpy as np\n",
    "import pandas\n",
    "from pymbar.timeseries import statisticalInefficiency"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cell sets up nglview with GUI controls (under development). When clicking on an atom, one can find all its details under the tab `Extra` and then `Picked`. We will often need the atom index in this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traj = mdtraj.load('traj.dcd', top='init.pdb')\n",
    "df = pandas.read_csv(\"scalars.csv\")\n",
    "view = nglview.show_mdtraj(traj, gui=True)\n",
    "view.add_representation(\"ball+stick\", selection=\"protein\")\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Inter-atomic distances\n",
    "\n",
    "The following example just computes two bond lengths as function of time, namely for the `N-CA` and `CA-C` bonds in the leucine residue. Note that all indices in Python start counting from zero, so the bond between the first two atoms in the PDB and DCD files is denoted by `[0, 1]`.\n",
    "\n",
    "The documentation of `compute_distances` can be found here: http://mdtraj.org/1.9.3/api/generated/mdtraj.compute_distances.html#mdtraj.compute_distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "distances = mdtraj.compute_distances(traj, [[0, 1], [1, 2]])\n",
    "print(distances.shape)\n",
    "print(distances[400, 0])\n",
    "plt.close(1)  # This is needed to rerun the code cell correctly\n",
    "fig, ax = plt.subplots(num=1)\n",
    "for counter, col in enumerate(distances.T):\n",
    "    ax.plot(df[\"Time (ps)\"], col, \".\", label=str(counter))\n",
    "ax.set_xlabel(\"Time [ps]\")\n",
    "ax.set_ylabel(\"Distances [nm]\")\n",
    "ax.legend(loc=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:#A03;font-size:14pt\">\n",
    "&#x270B; HANDS-ON! &#x1F528;\n",
    "</span>**\n",
    "\n",
    "> For the simple bond lengths shown above, the initialization (warming up of the protein during the first 5 picoseconds) seems to have little effect.\n",
    "> Try to find the atomic indices of strongly electrostatically interacting non-covalent atom pairs, e.g. the nitrogen in the lysine side changes interact with carboxilic groups in a glutamine and the c-terminus. Plot the distances between two pairs of electrostatically interacting atoms. Can you recognize the thermalization in these inter-atomic distances? Add two more distances, one for a hydrogen bond in an alpha helix and one outside an alpha helix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Valence and dihedral angles\n",
    "\n",
    "A similar analysis for the C-S-C angle, the first backbone dihedral angle and a side-chain dihedral angle, is carried out with [compute_angles](http://mdtraj.org/1.9.3/api/generated/mdtraj.compute_angles.html) and [compute_dihedrals](http://mdtraj.org/1.9.3/api/generated/mdtraj.compute_dihedrals.html) in the following code cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "angles = mdtraj.compute_angles(traj, [[171, 172, 173]])\n",
    "dihedrals = mdtraj.compute_dihedrals(traj, [[0, 1, 2, 21], [33, 36, 37, 39]])\n",
    "print(angles.shape)\n",
    "plt.close(2)  # This is needed to rerun the code cell correctly\n",
    "fig, ax = plt.subplots(num=2)\n",
    "ax.plot(df[\"Time (ps)\"], angles[:, 0]/np.pi*180, \".\", label='C-S-C angle')\n",
    "ax.plot(df[\"Time (ps)\"], dihedrals[:, 0]/np.pi*180,\n",
    "        \".\", label='N-CA-C-N angle (psi LEU1)')\n",
    "ax.plot(df[\"Time (ps)\"], dihedrals[:, 1]/np.pi*180,\n",
    "        \".\", label='CA-CB-CG-OD2 angle (chi2 ASP3)')\n",
    "ax.set_xlabel(\"Time [ps]\")\n",
    "ax.set_ylabel(\"Angles [deg]\")\n",
    "ax.legend(loc=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Valence angles in proteins are usually rather stiff while the ($\\psi$ and $\\phi$) dihedrals in the backbone and the $\\chi_n$ dihedrals in the side chains explain most of the conformational changes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Averages, fluctuations, error on the average\n",
    "\n",
    "Computing averages and fluctuations (standard deviations) of a time-dependent data series is easily carried out with Numpy. The following code cell contains some examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"average distance 0 [nm]\", distances[:, 0].mean())\n",
    "print(\"average distance 0 [nm]\", np.mean(distances[:, 0]))\n",
    "print(\"all average distances [nm]\", distances.mean(axis=0))\n",
    "print(\"all average distances [nm]\", np.mean(distances, axis=0))\n",
    "print(\"st.dev. distance 0 [nm]\", distances[:, 0].std())\n",
    "print(\"st.dev. distance 0 [nm]\", np.std(distances[:, 0]))\n",
    "print(\"all st.dev. distances [nm]\", distances.std(axis=0))\n",
    "print(\"all st.dev. distances [nm]\", np.std(distances, axis=0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantities from an MD simulation, of which an average is computed, are stochastic. Hence, any average over a finite number of steps is also a stochastic quantity, subject to an uncertainty. When computing an average over $N$ **uncorrelated** data points, the error on the average is\n",
    "\n",
    "$$\\sigma_{\\langle a \\rangle} = \\sqrt{\\frac{1}{N-1}\\sum_{i=1}^N (a_i - \\langle a \\rangle)^2}$$\n",
    "\n",
    "This formula is rarely useful for computing uncertainties on averages (over time) of quantities from MD simulations, because the values for subsequent time steps are often correlated. (See e.g. the inter-atomic distances for electrostatically interacting pairs of atoms.)\n",
    "\n",
    "One should correct the factor $1/(N-1)$ and replace $N$ by the number of uncorrelated data points in a time series. The function `num_independent` below estimates this number of independent samples with the block-averaging method. The implementation below is based on the description of the block average method in the book \"Computer Simulation of Liquids\" by M.P. Allen and D.J. Tildesley, section 6.4.1 on page 192."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def averror(values, num=None):\n",
    "    \"\"\"Compute the error on the average of (un)correlated samples.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    values\n",
    "        samples for which the error on the average is computed.\n",
    "    num\n",
    "        The number of independent samples. When not given, the\n",
    "        default is len(values), which would be correct for\n",
    "        uncorrelated samples.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    averror\n",
    "        The error on the averge.\n",
    "\n",
    "    \"\"\"\n",
    "    if num is None:\n",
    "        num = len(values)\n",
    "    elif num <= 1:\n",
    "        return np.nan\n",
    "    return np.std(values) / np.sqrt(num - 1)\n",
    "\n",
    "\n",
    "def num_independent(values, fignum=None):\n",
    "    \"\"\"Estimate the number of independent samples in a time series.\n",
    "\n",
    "    This implementation is based on the description in the book\n",
    "    of Allen and Tildesley. In addition to the book, this code\n",
    "    attempts to safeguard against balistic motion artifacts when\n",
    "    performing the extrapolation towards infinite block sizes.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    values\n",
    "        A time-correlated series, must be a numpy array with shape (N,).\n",
    "    fignum\n",
    "        A matplotlib figure number for plotting the inefficiency and the\n",
    "        model used for extrapolation.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    indep\n",
    "        The number of independent samples. When a constant input is provided,\n",
    "        or in a few pathological cases, the result will be 1.\n",
    "    quality\n",
    "        An empirical quality indicator to judge if the series is sufficiently\n",
    "        long for a reliable estimate of the number of independent samples.\n",
    "        This should be at least 5, preferrably more. When less than 5,\n",
    "        the time series should be made at least 2**(5-quality) times longer.\n",
    "\n",
    "    \"\"\"\n",
    "    if values.ndim != 1:\n",
    "        raise TypeError(\"Only one-dimensional arrays are supported.\")\n",
    "\n",
    "    # At least 16 blocks, then halve size at each iteration.\n",
    "    bss = []\n",
    "    ineffs = []\n",
    "    bs = len(values) // 16\n",
    "    while bs >= 1:\n",
    "        nb = len(values) // bs\n",
    "        blocks = values[:nb * bs].reshape(nb, bs)\n",
    "        avvar = blocks.mean(axis=1).var(ddof=1) / (nb - 1)\n",
    "        bss.insert(0, bs)\n",
    "        ineffs.insert(0, avvar)\n",
    "        bs //= 2\n",
    "\n",
    "    # Return in case of no blocks\n",
    "    if len(bss) < 1:\n",
    "        return 1, 0\n",
    "\n",
    "    # Convert to inefficiencies.\n",
    "    bss = np.array(bss)\n",
    "    ineffs = np.array(ineffs)\n",
    "    ineffs /= ineffs[0]\n",
    "\n",
    "    # Fit a simple model:  y = a * x / (a + x - 1)\n",
    "    # Discard decreasing part of the inefficiency, e.g. due to balistic motion.\n",
    "    i0 = ineffs.argmin()\n",
    "    y = ineffs[i0:] / ineffs[i0]\n",
    "    x = bss[i0:] / bss[i0]\n",
    "    # At least two data points are needed for the fit.\n",
    "    if len(y) >= 2:\n",
    "        dm = (y - x)/x\n",
    "        ev = (y * (1 - x))/x\n",
    "        a = np.dot(ev, dm) / np.dot(dm, dm)\n",
    "        ineff_limit = a * ineffs[i0]\n",
    "        optbs = a * bss[i0]\n",
    "        ineffs_model = ineff_limit * x / (a + x - 1)\n",
    "    else:\n",
    "        ineff_limit = 0\n",
    "        optbs = len(values)\n",
    "        ineffs_model = None\n",
    "\n",
    "    if ineff_limit <= 0 or ineff_limit > len(values):\n",
    "        num_indep = 1\n",
    "    else:\n",
    "        num_indep = len(values) / ineff_limit\n",
    "\n",
    "    # The quality measure is the minimum of the number of points\n",
    "    # used in the fit and the number of points after the optimal\n",
    "    # block size.\n",
    "    quality = min(len(y), (bss > optbs).sum())\n",
    "\n",
    "    if fignum is not None:\n",
    "        plt.close(fignum)\n",
    "        fig, ax = plt.subplots(num=fignum)\n",
    "        ax.set_title(\n",
    "            f\"Statistical inefficiency: {ineff_limit:.1f} ||\"\n",
    "            f\"Independent samples: {num_indep:.1f} ||\"\n",
    "            f\"Quality: {quality:d}\"\n",
    "        )\n",
    "        ax.plot(bss, ineffs, 'o', color='#aaaaaa')\n",
    "        ax.plot(bss[i0:], ineffs[i0:], 'ko')\n",
    "        if ineffs_model is not None:\n",
    "            ax.plot(bss[i0:], ineffs_model, '-', color='C0')\n",
    "            ax.axhline(ineff_limit, color='C0', ls=':')\n",
    "            ax.axvline(optbs, color='C0', ls=':')\n",
    "        ax.set_xscale(\"log\")\n",
    "        ax.set_xlabel(\"Block size\")\n",
    "        ax.set_ylabel(\"Statistical inefficiency\")\n",
    "\n",
    "    return num_indep, quality\n",
    "\n",
    "\n",
    "values = distances[:, 0]\n",
    "#values = angles[:, 0]\n",
    "#values = dihedrals[:, 0]\n",
    "num_indep, quality = num_independent(values, fignum=3)\n",
    "print(\"Quality:\", quality)\n",
    "print(\"Number of values:\", len(values))\n",
    "print(\"Number of independent values:\", num_indep)\n",
    "print(\"Error on the average:\", averror(values, num_indep))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:#A03;font-size:14pt\">\n",
    "&#x270B; HANDS-ON! &#x1F528;\n",
    "</span>**\n",
    "\n",
    "> Try estimating the error on the average of all internal coordinates computed so far. Does the number of independent samples correlate with the time dependence plotted previously?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Estimating the initialization phase (or burn-in)\n",
    "\n",
    "The error estimation in the previous section only works well when the provided time series is taken from an simulation in equilibrium, i.e. without initialization phase. One could try to estimate the length of the initialization or burn-in phase by visually inspecting the time-dependency in a plot. However, this is tedious, subjective and error-prone practice. In this section, we will employ the scheme proposed in [10.1021/acs.jctc.5b00784](https://doi.org/10.1021/acs.jctc.5b00784) to determine the initialization phase automatically. This scheme prescribes that the optimal initialization phase maximizes the number of independent samples in the remainder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scan_initialization(values):\n",
    "    \"\"\"Determine initialization phase.\n",
    "\n",
    "    This is a very approximate maximization of the number of independent\n",
    "    samples. A more fine-grained maximization is probably not that useful.\n",
    "    \"\"\"\n",
    "    ninit = 0\n",
    "    num_indep = num_independent(values)[0]\n",
    "    ninit_new = 1\n",
    "    while ninit_new < len(values):\n",
    "        num_indep_new = num_independent(values[ninit_new:])[0]\n",
    "        if num_indep_new > num_indep:\n",
    "            ninit = ninit_new\n",
    "            num_indep = num_indep_new\n",
    "        ninit_new = max(ninit_new + 1, int(ninit_new*1.1))\n",
    "    return ninit\n",
    "\n",
    "\n",
    "values = distances[:, 0]\n",
    "#values = angles[:, 0]\n",
    "#values = dihedrals[:, 0]\n",
    "#values = df[\"Potential Energy (kJ/mole)\"].to_numpy()\n",
    "ninit = scan_initialization(values)\n",
    "\n",
    "# Plot the time dependence and the initialization part\n",
    "plt.close(4)  # This is needed to rerun the code cell correctly\n",
    "fig, ax = plt.subplots(num=4)\n",
    "ax.plot(df[\"Time (ps)\"][:ninit], values[:ninit], \".\")\n",
    "ax.plot(df[\"Time (ps)\"][ninit:], values[ninit:], \".\")\n",
    "ax.axvline(df[\"Time (ps)\"][ninit], color='k')\n",
    "ax.set_xlabel(\"Time [ps]\")\n",
    "\n",
    "# Exhaustive plot of the number of independent samples as function\n",
    "# of the initialization time.\n",
    "tmp = np.array([num_independent(values[ninit:])\n",
    "               for ninit in range(len(values))])\n",
    "num_indeps, qualities = tmp.T\n",
    "plt.close(5)  # This is needed to rerun the code cell correctly\n",
    "fig, ax1 = plt.subplots(num=5)\n",
    "ax1 = plt.gca()\n",
    "ax1.set_xlabel('Initialization phase [1]')\n",
    "ax1.set_ylabel('Quality of the estimate', color=\"C0\")\n",
    "ax1.plot(qualities, color=\"C0\")\n",
    "ax1.axvline(ninit, color='k')\n",
    "ax1.tick_params(axis='y', labelcolor=\"C0\")\n",
    "ax2 = ax1.twinx()\n",
    "ax2.set_ylabel('Number of independent samples', color=\"C2\")\n",
    "ax2.plot(num_indeps, color=\"C2\", lw=1)\n",
    "ax2.tick_params(axis='y', labelcolor=\"C2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<span style=\"color:#A03;font-size:14pt\">\n",
    "&#x270B; HANDS-ON! &#x1F528;\n",
    "</span>**\n",
    "\n",
    "> Repeat the above code a few different time series and answer the following questions:\n",
    ">\n",
    "> - Visual inspection of the initialization phase can be too optimistic or pessimistic?\n",
    "> - What is the accuracy of an error estimate?\n",
    ">\n",
    "> An important consideration with the last question is that one has no information of unseen states in an MD simulation. For example, a finite MD simulation may be sampling a metastable conformation without relaxing to a more stable one. Any averages and error estimates derived from such a run are biased because they are not considering the thermodynamically dominant conformation.\n",
    "\n",
    "Note that in the reference [10.1021/acs.jctc.5b00784](https://doi.org/10.1021/acs.jctc.5b00784), an alternative (more robust but not as easy to understand) function is recommended to compute the number of independent samples in a time series. The original implementation of the authors can be found in [pymbar](https://pymbar.readthedocs.io/en/master/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(statisticalInefficiency(distances[:, 0]))\n",
    "print(statisticalInefficiency(angles[:, 0]))\n",
    "print(statisticalInefficiency(dihedrals[:, 0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Average of a dihedral angle\n",
    "\n",
    "The average of a dihedral angle can become meaningless because these angles are subject to a periodic boundary, i.e. $\\phi=35^\\circ$ and $\\phi=395^\\circ$ represent the same dihedral angles. As a consequence, dihedral angles exceeding $180^\\circ$ will continue at $-180^\\circ$. Such discontinuous jumps result in meaningless averages.\n",
    "\n",
    "This problem can be surmounted by computing a slightly different type of average, i.e. the average cosine and sine of the angle, followed by a conversion back into an angle. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.cos(dihedrals[:, 1]).mean()\n",
    "y = np.sin(dihedrals[:, 1]).mean()\n",
    "avdihed1 = np.arctan2(y, x)\n",
    "avdihed1_wrong = dihedrals[:, 1].mean()\n",
    "print(avdihed1/np.pi*180)\n",
    "print(avdihed1_wrong/np.pi*180)\n",
    "plt.close(6)  # This is needed to rerun the code cell correctly\n",
    "fig, ax = plt.subplots(num=6)\n",
    "ax.plot(\n",
    "    df[\"Time (ps)\"], dihedrals[:, 1]/np.pi*180,\n",
    "    \".\", label='CA-CB-CG-OD2 angle (chi2 ASP3)'\n",
    ")\n",
    "ax.axhline(avdihed1/np.pi*180, color='k')\n",
    "ax.axhline(avdihed1_wrong/np.pi*180, color='r')\n",
    "ax.set_xlabel(\"Time [ps]\")\n",
    "ax.set_ylabel(\"Angles [deg]\")\n",
    "ax.legend(loc=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the difference may be small but is certainly noticeable. The naive calculation of the average includes some samples close to $-180^\\circ$, whereas the alternative approach does not suffer from this issue. To understand why this works, consider a plot of the sine versus the cosine:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close(7)  # This is needed to rerun the code cell correctly\n",
    "fig, ax = plt.subplots(num=7)\n",
    "ax.plot(np.cos(dihedrals[:, 1]), np.sin(dihedrals[:, 1]), 'k+', alpha=0.1)\n",
    "ax.plot([x], [y], 'ro')\n",
    "ax.plot([0, x], [0, y], 'r-')\n",
    "ax.add_patch(mpatches.Arc(\n",
    "    [0, 0], 1, 1, angle=0, theta1=0, theta2=avdihed1/np.pi*180, color='r'\n",
    "))\n",
    "ax.axhline(0, color='k')\n",
    "ax.axvline(0, color='k')\n",
    "ax.set_xlim(-1.1, 1.1)\n",
    "ax.set_ylim(-1.1, 1.1)\n",
    "ax.set_aspect('equal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When plotting the dihedral angles as points on a circle, i.e. $x=\\cos\\phi$ and $y=\\sin\\phi$, there is no longer a discontinuity and the average is well-behaved. The angle is then derived from the average (red dot) by computing the angle with the X-axis, using the `arctan2` function. A few remarks:\n",
    "\n",
    "- The function `arctan2` takes the Y-coordinate as first argument.\n",
    "\n",
    "- This approach will not work when the dihedral angles uniformly distributed over the interval $[-180^\\circ,180^\\circ]$. In this case, the average cosine and sine are (nearly) zero and the angle with the X-axis is not (or ill) defined."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Histogram and free energy\n",
    "\n",
    "Histograms may also be convenient to characterize an internal coordinate. We will apply it here to the C-S-C angle, which is a relatively stiff mode. The angles are also [normally distributed](https://en.wikipedia.org/wiki/Normal_distribution), which can be shown by overlaying the histogram with the normal probability density.\n",
    "\n",
    "For reference, the force-field parameters for this angle can be found here:\n",
    "\n",
    "https://github.com/openmm/openmm/blob/master/wrappers/python/openmm/app/data/amber14/protein.ff14SB.xml#L3077\n",
    "\n",
    "and reads\n",
    "\n",
    "```\n",
    "<Angle angle=\"1.726130630222392\" k=\"518.816\" type1=\"protein-CT\" type2=\"protein-S\" type3=\"protein-CT\"/>\n",
    "```\n",
    "\n",
    "The equilibrium angle in degrees is $1.726130630222392 \\text{ rad}$, which corresponds to $98.9^\\circ$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.close(8)  # This is needed to rerun the code cell correctly\n",
    "fig, ax = plt.subplots(num=8)\n",
    "a = angles/np.pi*180  # for convenience\n",
    "bins = np.arange(85, 115, 1)\n",
    "ax.hist(a, bins, density=True)\n",
    "ava = a.mean()\n",
    "avs = a.std()\n",
    "x = np.linspace(80, 120, 79)\n",
    "# The following equation is the probability density of a\n",
    "# univariate normal distribution.\n",
    "y = np.exp(-(x - ava)**2 / (2*avs**2)) / np.sqrt(2*avs**2*np.pi)\n",
    "ax.plot(x, y)\n",
    "ax.axvline(ava, color='k')\n",
    "ax.set_title(f\"Average = {ava:.1f}, Std.Dev. = {avs:.1f}\")\n",
    "ax.set_xlabel(\"C-S-C angle [deg]\")\n",
    "ax.set_ylabel(\"Probability density [1/deg]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The average is slightly different from the equilibrium value, which may be due to other force-field terms influencing the local geometry.\n",
    "\n",
    "**<span style=\"color:#A03;font-size:14pt\">\n",
    "&#x270B; HANDS-ON! &#x1F528;\n",
    "</span>**\n",
    "\n",
    "> How can you judge if the average and the equilibrium are significantly different?\n",
    "\n",
    "The empirical probabilities from the histogram can be translated into free energies (up to a constant) by employing the following relationship:\n",
    "\n",
    "$$F(\\theta) = -k_B T \\log(p(\\theta))$$\n",
    "\n",
    "This relation can also be applied to the probability density to obtain a quadratic approximation of the free energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boltzmann = (1e-3 * 1.380649e-23 * 6.02214076e23)  # in kJ/mol\n",
    "temperature = 300\n",
    "# Force field parameters\n",
    "theta0_par = 1.726130630222392/np.pi*180  # rest angle converted to degrees\n",
    "# force constant converted to kJ/mol/degree^2\n",
    "fc_par = 518.816 / (180/np.pi)**2\n",
    "f0 = 0.5 * fc_par * (x - theta0_par)**2\n",
    "# Quadratic approximation\n",
    "fc = boltzmann * temperature / avs**2  # force constant\n",
    "f = 0.5 * fc * (x - ava)**2\n",
    "# Empiricial model, using histogram\n",
    "pe = np.histogram(a, bins, density=True)[0]\n",
    "with np.errstate(divide=\"ignore\"):\n",
    "    fe = -boltzmann * temperature * np.log(pe)\n",
    "fe -= fe.min()\n",
    "\n",
    "plt.close(9)  # This is needed to rerun the code cell correctly\n",
    "fig, ax = plt.subplots(num=9)\n",
    "ax.plot(x, f0, label='Force field valence angle term')\n",
    "ax.plot(x, f, label='Quadratic free energy model')\n",
    "ax.plot((bins[1:] + bins[:-1])/2, fe, label='Empirical free energy')\n",
    "ax.set_xlabel(\"C-S-C angle [deg]\")\n",
    "ax.set_ylabel(\"Free energy [kJ/mol]\")\n",
    "ax.legend(loc=0)\n",
    "\n",
    "# Finally, the force constant is printed in kJ/mol/rad^2, for comparison with the force-field parameter.\n",
    "print(fc * (180/np.pi)**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also here, the free-energy force constant deviates slightly from the AMBER force-field parameter. Again, this could be due to an influence from the environment of the C-S-C angle.\n",
    "\n",
    "**<span style=\"color:#A03;font-size:14pt\">\n",
    "&#x270B; HANDS-ON! &#x1F528;\n",
    "</span>**\n",
    "\n",
    "> Also estimate the error on the force constant."
   ]
  }
 ],
 "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"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}