{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "provenance": [] }, "kernelspec": { "name": "python3", "display_name": "Python 3" }, "language_info": { "name": "python" } }, "cells": [ { "cell_type": "markdown", "source": [ "#**Química Computacional - Aula Prática 13**\n", "\n", "---\n", "Bibliografia de Suporte:\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": [ "**Mecânica Estatística e Monte Carlo**\n", "\n", "\\\\\n", "\n", "**Exercício 1**\n", "\n", "Considere a distribuição de Boltzmann:\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) Mostre que a probabilidade relativa da ocorrência de dois estados do sistema com\n", "energias $E_0$ e $E_1$ pode ser calculada como:\n", "\n", "\\begin{equation}\n", " \\frac{P_1}{P_0} = e^{(-\\beta (E_1 - E_0)) }\n", "\\end{equation}\n", "\n", "b) Assumindo que o sistema se encontra em equilíbrio a uma temperatura $T$, mostre que\n", "a temperatura pode ser calculada como:\n", "\n", "\\begin{equation}\n", " T = \\frac{E_1 - E_0}{k_b ln(P_1/P_0)}\n", "\\end{equation}\n", "\n", "**Exercício 2**\n", "\n", "Considera a simulação através do método de Monte Carlo do modelo TIP-3P da água para um sistema com 256 moléculas à temperatura de 298.15 K. A energia potencial de uma determinada configuração, configuração “old”, é $U_o = -1.6918\\times 10^{-20} J$ . Após o deslocamento de uma molécula de água escolhida aleatoriamente obtém-se uma nova configuração (configuração “new”) com energia potencial $U_n = -1.3628\\times 10^{-20} J$\n", "\n", "\\\\\n", "a) Para decidir sobre a aceitação da nova configuração é gerado um número aleatório $η = 0.7501$ . Mostre, usando o algoritmo de Metropolis, que a nova configuração é rejeitada\n", "\n", "\\\\\n", "b) Após a rejeição da configuração “new” anterior gerou-se uma nova configuração, através do deslocamento de outra molécula de água escolhida aleatoriamente, com energia, $U_n = -1.6900\\times 10^{-20} J$ e um número aleatório $\\eta = 0.9503$ . Mostra que a nova configuração é aceite de acordo com o algoritmo de Metropolis.\n", "\n", "\\\\\n", "Nota 1:Se $U_n(\\vec{r}^{N'}) \\leq U_0(\\vec{r}^N)$, então $e^{-\\beta [U_n(\\vec{r}^{N'})-U_o(\\vec{r}^N)]} \\geq 1$: a configuração é aceite (probabilidade 1). Se $U_n(\\vec{r}^{N'}) >U_0(\\vec{r}^N)$ então a configuração é aceite com probabilidade $\\eta < e^{-\\beta [U_n(\\vec{r}^{N'})-U_o(\\vec{r}^N)]}$\n", "\n", "\\\\\n", "Nota 2: pode usar uma caixa de código para calcular os factor de Boltzmann $e^{-\\beta [U_n(\\vec{r}^{N'})-U_o(\\vec{r}^N)]}$\n", "\n", "\\\\\n", "Nota 3: Use $k_B = 1.38065 \\times 10^{-23} JK^{-1}$" ], "metadata": { "id": "CXoD_3r1u7QB" } }, { "cell_type": "code", "source": [ "import numpy as np\n", "\n" ], "metadata": { "id": "1tlBtM1mQ5m8" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "\n", "**Exercício 3**\n", "\n", "Construa agora um programa que decidida sobre a aceitação de novas configurações.\n", "Considere que $U_o = -1.6918\\times 10^{-20} J$ e foram geradas as seguintes configurações.\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", "Nota 1: Para cada nova configuração, gere um número aleatório e decida se a mesma é aceite ou rejeitada. Se a configuração for aceite, não se esqueça de a guardar para a próxima avaliação\n", "\n", "Nota 2: para gerar um número aleatório entre $[0,1]$ pode usar `np.random.rand()`\n" ], "metadata": { "id": "huip8S4hX5sb" } }, { "cell_type": "code", "source": [ "# Coloque o código aqui\n" ], "metadata": { "id": "ahIS8Lm_X19z" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercício 4**\n", "\n", "Vamos agora gerar um ensemble de 10000 configurações com a energia uniformemente distribuída de forma random no intervalo $[-1.9\\times 10^{-20} J, 0\\times 10^{-20} J]$. Construa um programa que decidida sobre a aceitação dessa novas configurações. Guarde as energias aceites para fazer um gráfico da sua distribuição. Calcule também a taxa de aceitação (número de configurações aceites/ numero de configurações totais) em %.\n", "\n", "\\\\\n", "Nota 1: para gerar $N$ números aleatórios entre $[a,b]$ pode usar `np.random.uniform(low=a, high=b, size=N)`\n", "\n", "\\\\\n", "Nota 2: Para armazena as energias aceites pode criar um vector `accepted_energies = [U_o]` que ao longo da iteração armazena os valores aceites\n", "\n", "```\n", "if i < boltzmann_factor:\n", " U_o = conf # aceitamos a configuração\n", " accepted_energies.append(U_o)\n", "```\n", "\n", "Podemos depois plotar um histograma\n", "\n", "\n", "```\n", "import matplotlib.pyplot as plt\n", "\n", "\n", "# plotar um histograma\n", "plt.hist([u * 1e20 for u in accepted_energies], bins=50)\n", "plt.xlabel(\"Energia ($10^{-20}$ J)\")\n", "plt.ylabel(\"Frequência\")\n", "plt.title(\"Distribuição das Energias Aceites\")\n", "plt.show()\n", "```" ], "metadata": { "id": "6g7DX2r8gm3s" } }, { "cell_type": "code", "source": [ "# Coloque o código aqui\n" ], "metadata": { "id": "q_lxDqZrgj41" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercício 5**\n", "\n", "Repita o exercício com as temperaturas $T = 500 K$ e $T = 90 K$. Compare a distribuição dos valores aceites e comente a taxa de aceitação." ], "metadata": { "id": "2AyFNhHxmHPA" } }, { "cell_type": "code", "source": [ "# Coloque o código aqui\n" ], "metadata": { "id": "XEfLVjWNXlgl" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercício 6**\n", "\n", "Considere agora um potencial harmonico usado para descrever uma ligação entre dois átomos\n", "\n", "$$\n", "U(x) = \\frac{1}{2}k(r-r_o)^2\n", "$$\n", "\n", "Os parâmetros do campo de forças GAFF (General Amber Force Field) para a ligação C(sp3)-C(sp3) são os seguintes\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 o Método de Monte Carlo/ Critério Metrópolis para calcular os valores médios da energia (e respectivos valor de $r$) para ligação C-C a 298.15 K. Indique se o valor médio da distância era o que esperava.\n", "\n", "Considere:\n", "\n", "```\n", "kB = 1.380649e-23 # J/K\n", "N = 100000 # número de configurações\n", "x_min = 0e-10 # limite inferior para x (em metros)\n", "x_max = 5e-10 # limite superior para x (em metros)\n", "```" ], "metadata": { "id": "BDp1UuFOoIdo" } }, { "cell_type": "code", "source": [ "# Coloque o código aqui\n", "\n" ], "metadata": { "id": "tEUD5dNu03MB" }, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "source": [ "**Exercício 7**\n", "\n", "Repita o exercício para $T = 90 K$ e $T = 500 K$. Compare as respectivas distribuições particularmente a distribuição das distâncias de ligação. Comente também o valor médio obtido a cada temperatura. Com base nos resultados anteriores, indique se à temperatura ambiente (~ 300K) poderia ser razoável considerarmos uma distância de ligação fixa durante uma simulação.\n" ], "metadata": { "id": "S6trYij-lOdE" } }, { "cell_type": "code", "source": [ "# Coloque o código aqui\n" ], "metadata": { "id": "H9wsvj3Klv7U" }, "execution_count": null, "outputs": [] } ] }