{
  "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 12**\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, Simulation and Molecular Modeling**\n",
        "\n",
        "\\\\\n",
        "\n",
        "**Exercise 1**\n",
        "\n",
        "Applying the principle of equipartition of kinetic energy to the instantaneous kinetic energy of the system, it is possible to obtain a desired instantaneous temperature $T_d$, scaling the components of the velocities of the $N$ particles of the system\n",
        "\n",
        "\\begin{equation}\n",
        "v'_{i \\alpha} = v_{i \\alpha} \\sqrt{\\frac{T_d}{T_a}}\n",
        "\\end{equation}\n",
        "\n",
        "Show that the reduction of the components $α = x, y, z$ of the velocity of a system of $N$ monatomic particles through the previous equation actually produces the desired instantaneous temperature ($T_d$).\n",
        "\n",
        "Suggestion: start by writing the expression for the kinetic energy as a function of the desired velocities.\n",
        "\n",
        "\\\\\n",
        "**Exercise 2**\n",
        "\n",
        "In the absence of external forces, linear momentum is conserved. Therefore, the initial velocities must be scaled so that the total linear momentum of the system is initially zero. This can be done using the equation:\n",
        "\n",
        "\\begin{equation}\n",
        "v'_{i \\alpha} (0) = v_{i \\alpha} (0) - \\frac{1}{N} \\sum^{N}_{j=1} v_{j \\alpha} (0) \\ \\text{with $\\alpha = x,y,x$}\n",
        "\\end{equation}\n",
        "\n",
        "Show that the $α$ component of the total linear momentum of a system of $N$ atoms with atomic mass $m$, in which the velocity components of each particle have been reduced using the previous equation, is in fact zero\n",
        "\n",
        "\\\\\n",
        "**Exercise 3**\n",
        "\n",
        "The Lennard-Jones potential\n",
        "\n",
        "\\\\\n",
        "\\begin{equation}\n",
        "u(r) = 4 \\epsilon \\left[ \\left(\\frac{\\sigma}{r}\\right)^{12} - \\left(\\frac{\\sigma}{r}\\right)^{6}\\right]\n",
        "\\end{equation}\n",
        "\n",
        "can be used to describe the Pauli repulsions and London attractions between particles.\n",
        "\n",
        "\\\n",
        "On the other hand, the magnitude of the force felt by particle $i$ due to the interaction with the other $N-1$ particles in the system is given by\n",
        "\n",
        "\\\\\n",
        "\\begin{equation*}\n",
        "F_{i} = - \\sum^{N}_{j \\neq i} \\frac{d \\ u(r_{ij})}{d \\ r_{ij}}\n",
        "\\end{equation*}\n",
        "\n",
        "Calculate the magnitude of the force felt by particle $i$ due to the interaction with the other $N-1$ under the effect of a Lennard-Jones potential\n",
        "\n",
        "Note 1: You will need the rule of exponents $\\frac{d}{dr} (r^{n}) = n r^{n-1}$\n",
        "\n",
        "\\\\\n",
        "**Exercise 4**\n",
        "\n",
        "Using the Verlet method, the calculation of the positions at time $(t+ Δ t)$ is given by\n",
        "\n",
        "$$\n",
        "\\mathbf{r}(t + \\Delta t) = 2\\mathbf{r}(t) - \\mathbf{r}(t - \\Delta t) + \\frac{1}{m} \\mathbf{F}(t)\\Delta t^2\n",
        "$$\n",
        "\n",
        "\\\\\n",
        "Calculate the $x$ component of the position vector of an atom $i$ at time $t_1 = ( t_0 + ∆t)$ of a molecular dynamics simulation, knowing that:\n",
        "\n",
        "\\\\\n",
        "$x_i(t_0) = 0$, $v_{ix}(t_0) = 1.0$, $a_{ix}(t_0) = 0.5$ and $Δt= 0.001$\n",
        "\n",
        "\\\\\n",
        "Recall that we can estimate $x(t-\\Delta t)$ when $t=t_0$ using the Taylor series truncated at the second-order term\n",
        "\n",
        "\\\\\n",
        "$$\n",
        "x(t_0-\\Delta t) = x(t_0) - v_x(t_0)\\Delta t\n",
        "$$\n",
        "\n",
        "\\\\\n",
        "Assume a reduced system of units (i.e., all quantities are dimensionless)"
      ],
      "metadata": {
        "id": "CXoD_3r1u7QB"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 5**\n",
        "\n",
        "Consider a particle moving under a constant force, for example, a free fall under the action of gravity\n",
        "\n",
        "$$\n",
        "𝐹(𝑥)=−mg\n",
        "$$\n",
        "\n",
        "Assume a reduced system of units (i.e., all quantities are dimensionless) with the following values:\n",
        "\n",
        "* Mass: $1$\n",
        "* Initial position: $x_0 = 10.0$\n",
        "* Initial velocity: $v_0 = 0.0$\n",
        "* Acceleration due to gravity: $g = a = 9.8$\n",
        "* Time-step: $Δt = 0.001$\n",
        "* Total simulation time: $20$\n",
        "\n",
        "Simulate the trajectory of the particle over time using the Verlet method. Plot the position over time.\n",
        "\n",
        "Note: Use the following code as a template\n",
        "\n",
        "```\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Set the parameters\n",
        "g = -9.8 # acceleration\n",
        "m = 1.0 # mass\n",
        "dt = 0.01 # time step\n",
        "t_tot=20 # total time\n",
        "\n",
        "# the number of points is given by t_tot/dt\n",
        "n_steps = int(t_tot / dt)\n",
        "\n",
        "# Initializing the arrays\n",
        "#\n",
        "# Let's initialize a time vector from 0 to 20 containing n_step points\n",
        "t = np.linspace(0, t_tot, n_steps)\n",
        "# let's initialize a position vector x with the same size\n",
        "# let's initialize a velocity vector v with the same size\n",
        "x = np.zeros(n_steps)\n",
        "v = np.zeros(n_steps)\n",
        "\n",
        "# Initial conditions:\n",
        "x[0] = 10.0 # initial position\n",
        "v[0] = 0.0 # initial velocity\n",
        "a = g # acceleration (constant)\n",
        "\n",
        "# In the Verlet method, we need to estimate x(t-Δt) when t=t_0\n",
        "# using an approximation: Taylor series of x(t−∆t) truncated at the second-order term):\n",
        "# x(t_0-∆t) = x(t_0) - v_x(t_0)Δt\n",
        "\n",
        "# Since x(t0) = x[0] (initial position)\n",
        "# x(t0−∆t) = x[-1] (position at previous step)\n",
        "# Estimate x[-1]\n",
        "x_minus_1 = x[0] - v[0]*dt\n",
        "\n",
        "# Verlet integration method\n",
        "# x(t+∆t) = 2 x(t) - x(t-∆t) + a ∆t^2\n",
        "#\n",
        "# Once we have assigned the initial values ​​x[0], v[0], we will iterate from 1 to n_steps\n",
        "# to calculate x[1], x[2], ..., x[n_steps], i.e.\n",
        "# we will calculate the positions at discrete points of the trajectory\n",
        "\n",
        "for i in range(1, n_steps):\n",
        "# calculate the position for i = t+∆t\n",
        "x[i] = 2*x[i-1] - x_minus_1 + (a * dt**2)\n",
        "# update the previous position, i.e.\n",
        "# change the value of the variable x_minus_1\n",
        "x_minus_1 = x[i-1]\n",
        "```\n"
      ],
      "metadata": {
        "id": "nPL7Z3LX0Bxu"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Set the parameters\n",
        "g = -9.8 # acceleration\n",
        "m = 1.0 # mass\n",
        "dt = 0.001 # time step\n",
        "t_tot=20 # total time\n",
        "\n",
        "# Put your code here\n"
      ],
      "metadata": {
        "id": "UGKCYlr_1LAV"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 6**\n",
        "\n",
        "Using this algorithm, we can calculate the velocities in order of time using:\n",
        "\n",
        "\\\\\n",
        "\\begin{equation}\n",
        "\\mathbf{v}(t) ≈ \\frac{\\mathbf{r}(t + \\Delta t) - \\mathbf{r}(t - \\Delta t)}{2\\Delta t}\n",
        "\\end{equation}\n",
        "\n",
        "Adapt the previous program to calculate the velocity and represent the value of the velocity as a function of time. Note that the velocities are calculated at time $t$ and not at $t + ∆t$"
      ],
      "metadata": {
        "id": "2lIbKn5N86db"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# Coloque o seu código aqui"
      ],
      "metadata": {
        "id": "IzJkK-n_9aCj"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 7**\n",
        "\n",
        "Now consider harmonic potential\n",
        "\n",
        "$$\n",
        "𝐹(𝑥)=−𝑘𝑥\n",
        "$$\n",
        "\n",
        "where $k$ is the force constant. Assume a reduced system of units (i.e., all quantities are dimensionless) with the following values:\n",
        "\n",
        "* Mass: $1$\n",
        "* Initial position: $x_0 = 1.0$\n",
        "* Initial velocity: $v_0 = 0.0$\n",
        "* Time step: $Δt = 0.001$\n",
        "* Total simulation time: $20$\n",
        "\n",
        "Simulate the trajectory over time using the Verlet method adapting the previous code. Also check the effect of the time step (e.g., repeat the procedure using $Δt = 0.5$)\n"
      ],
      "metadata": {
        "id": "rlQ4hWgYz_I8"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# coloque o seu código aqui"
      ],
      "metadata": {
        "id": "sxS2NX3UZhUm"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 8**\n",
        "\n",
        "As an alternative to Verlet's method, we can use the Leapfrog method to integrate the equations of motion:\n",
        "\n",
        "\\\\\n",
        "$$\n",
        "\\begin{equation}\n",
        "    \\begin{aligned}\n",
        "&\\mathbf{v}(t+{{\\frac{1}{2}}{{\\Delta t}}}) = \\mathbf{v}(t-{{\\frac{1}{2}}{{\\Delta t}}})+\\frac{{{\\Delta t}}}{m}\\mathbf{F}(t) \\\\\n",
        "&\\mathbf{r}(t+{{\\Delta t}}) = \\mathbf{r}(t)+{{\\Delta t}} \\ \\mathbf{v}(t+{{\\frac{1}{2}}{{\\Delta t}}})\n",
        "\\end{aligned}\n",
        "\\end{equation}\n",
        "$$\n",
        "\n",
        "\\\\\n",
        "Note that once again you can obtain an estimate of the velocity $v(\\frac{Δt}{2})$ when $t=0$ using a truncated Taylor series\n",
        "\n",
        "\\\\\n",
        "$$\n",
        "v(t+\\frac{Δt}{2}) = v(t)+ \\frac{1}{2}a(t)Δt\n",
        "$$\n",
        "\n",
        "\\\\\n",
        "Repeat exercise 5 using this integrator. Compare the trajectory with that obtained with the Verlet method"
      ],
      "metadata": {
        "id": "lx-q1CeYiZhx"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Coloque aqui o seu código\n"
      ],
      "metadata": {
        "id": "DmsDAW3xIEIq"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "**Exercise 9**\n",
        "\n",
        "Now use the Leapfrog method to simulate the trajectory of exercise 5. Check the effect of the time-step once more."
      ],
      "metadata": {
        "id": "S2vS-3tJOvhA"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Coloque aqui o seu código\n"
      ],
      "metadata": {
        "id": "a6Nr33ZNjO3r"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}