{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "provenance": [] }, "kernelspec": { "name": "python3", "display_name": "Python 3" }, "language_info": { "name": "python" } }, "cells": [ { "cell_type": "markdown", "source": [ "#**Computational Chemistry - Practical Class 13**\n", "\n", "---\n", "Supporting Bibliography:\n", "\n", "Mike P. Allen and Dominic J. Tildesley, ”Computer Simulation of Liquids”,\n", "Clarendon Press, Oxford, 1997\n", "\n", "---" ], "metadata": { "id": "8crRlPgWWMi6" } }, { "cell_type": "markdown", "source": [ "**Statistical Mechanics and Monte Carlo**\n", "\n", "\\\\\n", "\n", "**Exercise 1**\n", "\n", "Consider the Boltzmann distribution:\n", "\n", "\\begin{equation}\n", "P_i = \\frac{e^{-E_i \\beta}}{\\sum_{i} e^{E_i \\beta}} = \\frac{e^{-E_i \\beta}}{Q(N,V,T)}\n", "\\end{equation}\n", "\n", "a) Show that the relative probability of the occurrence of two states of the system with\n", "energies $E_0$ and $E_1$ can be calculated as:\n", "\n", "\\begin{equation}\n", "\\frac{P_1}{P_0} = e^{(-\\beta (E_1 - E_0)) }\n", "\\end{equation}\n", "\n", "b) Assuming that the system is in equilibrium at a temperature $T$, show that\n", "the temperature can be be calculated as:\n", "\n", "\\begin{equation}\n", " T = \\frac{E_1 - E_0}{k_b ln(P_1/P_0)}\n", "\\end{equation}\n", "\n", "**Exercise 2**\n", "\n", "Consider the simulation using the Monte Carlo method of the TIP-3P model of water for a system with 256 molecules at a temperature of 298.15 K. The potential energy of a certain configuration, the “old” configuration, is $U_o = -1.6918\\times 10^{-20} J$ . After the displacement of a randomly chosen water molecule, a new configuration is obtained (“new” configuration) with potential energy $U_n = -1.3628\\times 10^{-20} J$\n", "\n", "\\\\\n", "a) To decide whether to accept the new configuration, a random number $η = 0.7501$ is generated. Show, using the Metropolis algorithm, that the new configuration is rejected\n", "\n", "\\\\\n", "b) After rejecting the previous “new” configuration, a new configuration was generated by moving another randomly chosen water molecule with energy $U_n = -1.6900\\times 10^{-20} J$ and a random number $\\eta = 0.9503$ . Show that the new configuration is accepted according to the Metropolis algorithm.\n", "\n", "\\\\\n", "Note 1: If $U_n(\\vec{r}^{N'}) \\leq U_0(\\vec{r}^N)$, then $e^{-\\beta [U_n(\\vec{r}^{N'})-U_o(\\vec{r}^N)]} \\geq 1$: the configuration is accepted (probability 1). If $U_n(\\vec{r}^{N'}) >U_0(\\vec{r}^N)$ then the configuration is accepted with probability $\\eta < e^{-\\beta [U_n(\\vec{r}^{N'})-U_o(\\vec{r}^N)]}$\n", "\n", "\\\\\n", "Note 2: You can use a code box to calculate the Boltzmann factors $e^{-\\beta [U_n(\\vec{r}^{N'})-U_o(\\vec{r}^N)]}$\n", "\n", "\\\\\n", "Note 3: Use $k_B = 1.38065 \\times 10^{-23} JK^{-1}$" ], "metadata": { "id": "CXoD_3r1u7QB" } }, { "cell_type": "code", "source": [ "# Put the code here\n", "\n" ], "metadata": { "id": "1tlBtM1mQ5m8" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercise 3**\n", "\n", "Now build a program that decides whether to accept new configurations.\n", "\n", "Consider that $U_o = -1.6918\\times 10^{-20} J$ and the following configurations were generated.\n", "\n", "$U_{n1} = -1.6900\\times 10^{-20} J$\n", "\n", "$U_{n2} = -1.6950\\times 10^{-20} J$\n", "\n", "$U_{n3} = -1.3628\\times 10^{-20} J$\n", "\n", "$U_{n4} = -1.7000\\times 10^{-20} J$\n", "\n", "$U_{n4} = -0.2\\times 10^{-20} J$\n", "\n", "Note 1: For each new configuration, generate a random number and decide whether it is accepted or rejected. If the configuration is accepted, don't forget to save it for the next evaluation\n", "\n", "Note 2: To generate a random number between $[0,1]$ you can use `np.random.rand()`" ], "metadata": { "id": "huip8S4hX5sb" } }, { "cell_type": "code", "source": [ "# Put the code here\n" ], "metadata": { "id": "ahIS8Lm_X19z" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercise 4**\n", "\n", "Let us now generate an ensemble of 10000 configurations with energy uniformly distributed randomly in the interval $[-1.9\\times 10^{-20} J, 0\\times 10^{-20} J]$. Build a program that decides whether to accept these new configurations. Save the accepted energies to make a graph of their distribution. Also calculate the acceptance rate (number of accepted configurations/number of total configurations) in %.\n", "\n", "\\\\\n", "Note 1: to generate $N$ random numbers between $[a,b]$ you can use `np.random.uniform(low=a, high=b, size=N)`\n", "\n", "\\\\\n", "Note 2: To store the accepted energies you can create a vector `accepted_energies = [U_o]` that stores the accepted values ​​throughout the iteration\n", "\n", "```\n", "if i < boltzmann_factor:\n", "U_o = conf # we accept the configuration\n", "accepted_energies.append(U_o)\n", "```\n", "\n", "We can then plot a histogram\n", "\n", "```\n", "import matplotlib.pyplot as plt\n", "\n", "# plot a histogram\n", "plt.hist([u * 1e20 for u in accepted_energies], bins=50)\n", "plt.xlabel(\"Energy ($10^{-20}$ J)\")\n", "plt.ylabel(\"Frequency\")\n", "plt.title(\"Distribution of Accepted Energies\")\n", "plt.show()\n", "```" ], "metadata": { "id": "6g7DX2r8gm3s" } }, { "cell_type": "code", "source": [ "# Put the code here\n" ], "metadata": { "id": "q_lxDqZrgj41" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercise 5**\n", "\n", "Repeat the exercise with temperatures $T = 500 K$ and $T = 90 K$. Compare the distribution of accepted values ​​and comment on the acceptance rate." ], "metadata": { "id": "2AyFNhHxmHPA" } }, { "cell_type": "code", "source": [ "# Put the code here\n" ], "metadata": { "id": "XEfLVjWNXlgl" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercise 6**\n", "\n", "Now consider a harmonic potential used to describe a bond between two atoms\n", "\n", "$$\n", "U(x) = \\frac{1}{2}k(r-r_o)^2\n", "$$\n", "\n", "The parameters of the GAFF (General Amber Force Field) force field for the C(sp3)-C(sp3) bond are as follows\n", "\n", "$k= 300.9~kcal/mol/Å^2 ≈ 209.06 N/m$\n", "\n", "$r_0 =1.5375 Å = 1.5375 \\times 10^{−10} m$\n", "\n", "Use the Monte Carlo Method/Metropolis Criterion to calculate the mean energy values ​​(and respective $r$ values) for the C-C bond at 298.15 K. Indicate whether the mean distance value was what you expected.\n", "\n", "Consider:\n", "\n", "```\n", "kB = 1.380649e-23 # J/K\n", "N = 100000 # number of configurations\n", "x_min = 0e-10 # lower bound for x (in meters)\n", "x_max = 5e-10 # upper bound for x (in meters)\n", "```" ], "metadata": { "id": "BDp1UuFOoIdo" } }, { "cell_type": "code", "source": [ "# Put your code here\n", "\n" ], "metadata": { "id": "tEUD5dNu03MB" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercise 7**\n", "\n", "Repeat the exercise for $T = 90 K$ and $T = 500 K$. Compare the respective distributions, particularly the distribution of bond distances. Also comment on the mean value obtained at each temperature. Based on the previous results, indicate whether at room temperature (~ 300 K) it would be reasonable to consider a fixed bond distance during a simulation.\n" ], "metadata": { "id": "S6trYij-lOdE" } }, { "cell_type": "code", "source": [ "# Put your code here\n" ], "metadata": { "id": "H9wsvj3Klv7U" }, "execution_count": null, "outputs": [] } ] }