{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Molecular Dynamics" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "We have introduced the classical potential models, and have derived and showen some of their basic properties.\n", "Now we can use these potential models to look at the dynamics of the system. \n", "\n", "## Force and acceleration\n", "The particles that we study are classical in nature, therefore we can apply classical mechanics to rationalise their dynamic behaviour. \n", "For this the starting point is **Newton's second law of motion**, \n", "\n", "$$ \\mathbf{f} = m\\mathbf{a}, $$\n", "\n", "where $\\mathbf{f}$ is the force vector on an atom of mass, $m$, with an acceleration vector, $\\mathbf{a}$. \n", "The force, $f$, between two particles, at a position $r$, can be found from the interaction energy, $E(r)$, \n", "\n", "$$ f = \\dfrac{-\\partial E(r)}{\\partial r}.$$\n", "\n", "Which is to say that the force is the negative of the first derivative of the energy with respect to the postion of the particles.\n", "The Python code below creates a new function that is capable of calculating the force from the Lennard-Jones potential. \n", "The force on the atoms is then plotted." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEPCAYAAABsj5JaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xt4XPV95/H3d2Y0Gt0lS7Ity5ZlY2Ns7sbhUgIhOCTQXCApdCHJk3RDSzcbdp9uNm3JbpdNaZo2bTds2pC2PLmxbBpIKHSdFJob5NYGJ7a5+op8v1uWLMm6a6Tv/jFHlqxItjTSzBlpPq/n0TPnnPlJ58s8WB/9fr9zfsfcHRERkamKhF2AiIjMTgoQERFJiwJERETSogAREZG0KEBERCQtChAREUmLAkRERNKiABERkbQoQEREJC2xsAvIpJqaGm9sbAy7DBGRWWXz5s0n3b32fO3mdIA0NjayadOmsMsQEZlVzGz/ZNppCEtERNKiABERkbQoQEREJC0KEBERSYsCRERE0qIAERGRtGQ9QMzsVjPbaWZNZvbAOO8XmtmTwfsbzaxxzPsNZtZpZp/IVs0iIvKrshogZhYFHgFuA9YA95jZmjHN7gVOufsK4GHgs2Pe/xzwXCbrPNLWw+e+t5O9J7syeRoRkVkt2z2Qq4Emd9/j7v3AE8DtY9rcDjwWbD8FrDczAzCzO4C9wNZMFtna1c9fP9/EruOnM3kaEZFZLdsBUg8cHLV/KDg2bht3TwLtQLWZlQJ/CPxxpousLC4AoK27P9OnEhGZtWbTJPqngIfdvfNcjczsPjPbZGabmpub0zpRVXEcgLbugbS+X0QkH2R7LazDwJJR+4uDY+O1OWRmMaACaAGuAe40s78AKoEhM+t19y+M/mZ3fxR4FGDdunWeTpHF8SgFUeOUAkREZELZDpBfAivNbBmpoLgbeP+YNhuADwM/B+4Ennd3B24YbmBmnwI6x4bHTDEzKovjtPdoCEtEZCJZDRB3T5rZ/cB3gSjwFXffamYPAZvcfQPwZeBxM2sCWkmFTNZVFhVoCEtE5Byyvpy7uz8LPDvm2IOjtnuBu87zMz6VkeJGqSqOc0qT6CIiE5pNk+hZVVGsHoiIyLkoQCagISwRkXNTgEygqiROmybRRUQmpACZQEVRAb0DQ/QODIZdiohITlKATEA3E4qInJsCZALDy5noSiwRkfEpQCYwsh6WeiAiIuNRgEygsig1hKW70UVExqcAmUBVyfAQlnogIiLjUYBMYLgHoiEsEZHxKUAmUBSPUhiL6JkgIiITUICcQ6WWMxERmZAC5By0oKKIyMQUIOdQUVRAW496ICIi41GAnENqCEs9EBGR8ShAzqGqOK45EBGRCShAzqGiODWElXqiroiIjKYAOYeq4jj9ySF6tCKviMivUICcQ2WR1sMSEZmIAuQcKoMl3XUpr4jIr1KAnMPwirzt6oGIiPwKBcg5VJ3pgShARETGUoCcw5lngmhJdxGRX6EAOYcKTaKLiExIAXIOiYIoRQVR3Y0uIjIOBch5VGlFXhGRcSlAzqOiOK5JdBGRcShAzqOyqEDPRRcRGYcC5DyqSgrUAxERGYcC5DwqirQir4jIeBQg51EVPBNEK/KKiJxNAXIelcUFJIecrn6tyCsiMpoC5DzOLKjYpYl0EZHRsh4gZnarme00syYze2Cc9wvN7Mng/Y1m1hgcv9rMXg6+XjGz92aj3uEl3dv1bHQRkbNkNUDMLAo8AtwGrAHuMbM1Y5rdC5xy9xXAw8Bng+OvA+vc/QrgVuDvzSyW6ZqrSlI9EE2ki4icLds9kKuBJnff4+79wBPA7WPa3A48Fmw/Baw3M3P3bndPBscTQFZmtYd7IHomiIjI2bIdIPXAwVH7h4Jj47YJAqMdqAYws2vMbCvwGvAfRgXKGWZ2n5ltMrNNzc3N0y54XtADaensm/bPEhGZS2bVJLq7b3T3i4E3AZ80s8Q4bR5193Xuvq62tnba56wqjhONGCc71QMRERkt2wFyGFgyan9xcGzcNsEcRwXQMrqBu28HOoFLMlZpIBIx5pXEOakeiIjIWbIdIL8EVprZMjOLA3cDG8a02QB8ONi+E3je3T34nhiAmS0FLgL2ZaPomtJCBYiIyBgZv4ppNHdPmtn9wHeBKPAVd99qZg8Bm9x9A/Bl4HEzawJaSYUMwJuBB8xsABgC/qO7n8xG3TWlcZo1hCUicpasBgiAuz8LPDvm2IOjtnuBu8b5vseBxzNe4DhqSwvZ09wVxqlFRHLWrJpED0tNWWoIS+thiYiMUIBMQk1pnL7kEJ19v3LVsIhI3lKATEJNaSGALuUVERlFATIJIwGiK7FERIYpQCbhTICcVoCIiAxTgExCTVlqORP1QERERihAJmFecRwzdC+IiMgoCpBJiEUjzCvWciYiIqMpQCapprRQcyAiIqMoQCappkw9EBGR0RQgk5RaUFFzICIiwxQgk6QVeUVEzqYAmaSa0kK6+wfp7tdyJiIioACZtJrS4F6Q0xrGEhEBBcik1ZSl7kZv1jCWiAigAJm0Wq2HJSJyFgXIJGlBRRGRsylAJqlacyAiImdRgExSQTRCZXGBeiAiIgEFyBToXhARkREKkCmoKdVyJiIiwxQgU6DlTERERihApkAr8oqIjFCATEFtWSGn+5L0DgyGXYqISOgUIFNwZjkTzYOIiChApmLkZkLNg4iIKECm4EyAaB5EREQBMhXDCypqCEtERAEyJcMLKh7vUICIiChApiAei1BTWsjR9p6wSxERCZ0CZIoWVSY40t4bdhkiIqFTgEzRwvIEx9QDERHJfoCY2a1mttPMmszsgXHeLzSzJ4P3N5pZY3D8FjPbbGavBa83Z7t2gEWVRRxtUw9ERCSrAWJmUeAR4DZgDXCPma0Z0+xe4JS7rwAeBj4bHD8JvNvdLwU+DDyenarPVleR4HRfktO9A2GcXkQkZ0w5QMzsD6dxvquBJnff4+79wBPA7WPa3A48Fmw/Baw3M3P3l9z9SHB8K1BkZoXTqCUtdZVFABzVPIiI5LnY+RqY2TdH7wJXMNIrmKp64OCo/UPANRO1cfekmbUD1aR6IMN+A9ji7lm/nrauIgGkAuTCBWXZPr2ISM44b4AAHe7+28M7Zva3GaznvMzsYlIB9vYJ3r8PuA+goaFhxs9/JkDaNJEuIvltMkNYfzpm/79P43yHgSWj9hcHx8ZtY2YxoAJoCfYXA88AH3L33eOdwN0fdfd17r6utrZ2GqWOb0F5AjN0Ka+I5L3zBoi77wUws5pgv3Ua5/slsNLMlplZHLgb2DCmzQZSk+QAdwLPu7ubWSXwz8AD7v6v06hhWgqiEeaXFaoHIiJ5byqT6F+Z7sncPQncD3wX2A580923mtlDZvaeoNmXgWozawI+Dgxf6ns/sAJ40MxeDr7mT7emdCysKNIkuojkvcnMgQyzmTihuz8LPDvm2IOjtnuBu8b5vk8Dn56JGqZrUUWCXcdPh12GiEioptID8YxVMcvUBT0Qd30kIpK/phIgM9IDmQsWVSbo7h+koycZdikiIqGZSoB8MmNVzDILg0t5j2hNLBHJY5MOEHd/feyxYN2qu8zsn2a2rNxWV5G6G/2YJtJFJI+ls5SJmdl6M/sqcAR4J/C1mS4sly2qVA9ERGQyS5l8GPhN4EVSN/W9D3iZ1DpWN7v7b2WywFw0vyxBNGJalVdE8tpkLuP9OLAeOAH8BLjJ3Q8AmNmfZ7C2nBWNGPPLCtUDEZG8NpkA+RtSq+P+MVAKvGBmLwPfII8fSFVXkdAciIjktfMGiLt/CfjS8L6Z/QFwM/ABoMzMHgeecfenM1ZlDqqrLGLbkY6wyxARCc2kexBm9oiZXe8pP3T3j5Baev07wIcyVmGOWlSR4Ehbj24mFJG8NZUhqF3AX5nZPjP7CzO7wt173f1Jd78jUwXmqrqKIvqSQ5zq1pMJRSQ/TeU+kM+7+3XAW0gtr/5VM9thZv/TzFZmrMIcNfJgKU2ki0h+mvIkuLvvd/fPuvuVwD3AHcCOGa8sx515tK0u5RWRPJXOjYQxM3u3mX0deA7YSerekLyySD0QEclzk17O3cxuIdXj+HXgF6RuJLzP3bsyVFtOqyktpCBqHFYPRETy1FSeB/JJ4B+A/+rupzJUz6wRiRhL5hVzoDUv81NEZPIB4u43Z7KQ2aixuoR9J7vDLkNEJBTpLqb4QTN7MNhvMLOrZ7603Le0upj9LV26F0RE8lI6S5F8EbiO1HwIwGngkRmraBZprC6hq3+Qk539YZciIpJ16QTINe7+MaAXIJgPic9oVbPE0upiAPa3aB5ERPJPOgEyYGZRgmekm1ktMDSjVc0SjdUlAOxr0TyIiOSfdALkr4FngPlm9qfAz4DPzGhVs0R9VRHRiKkHIiJ5aSr3gVwHvOjuXzezzaSeEWLAHe6+PVMF5rKCaITFVUXqgYhIXprKfSAfAh4xs13AvwD/6O7HMlPW7LG0ukQ9EBHJS1MZwup397XAQ0AV8DUz+7mZfcbMbgzmRfJOY3Uxe0/qUl4RyT9TCZAbg9fH3f1hd7+V1IOlfgbcBWyc6eJmg6XVJZzuTdKmZd1FJM9MZQjrh2b2c2ChmX0EeAV43d2fBZ7NSHWzQGNwKe++li6qSvLyamYRyVNTeR7IJ4APAoPAMuB/AK+b2VYzezJD9eW8pcGlvPs1kS4ieWYqPRDcfbeZvc3ddw0fM7NS4JIZr2yWWDKvCLNUD0REJJ9MKUAARodHsN8JvDhjFc0yhbEoiyqK1AMRkbyTzo2EMkZjTbF6ICKSdxQgMyB1L4h6ICKSXxQgM6CxupjWrn7ae3Qpr4jkj6wHiJndamY7zazJzB4Y5/1CM3syeH+jmTUGx6vN7AUz6zSzL2S77nMZvhLrgHohIpJHshogwd3qjwC3AWuAe8xszZhm9wKn3H0F8DDw2eB4L6lLhz+RpXInbWRVXs2DiEj+yHYP5Gqgyd33uHs/8ARw+5g2twOPBdtPAevNzNy9y91/RvAcklzSME/PBRGR/JPtAKkHDo7aPxQcG7eNuyeBdqB6sicws/vMbJOZbWpubp5muZNTFI9SV5Fgd7MCRETyx5ybRHf3R919nbuvq62tzdp5Vy0sY8ex01k7n4hI2LIdIIeBJaP2FwfHxm1jZjGgAmjJSnXTsGphGbtPdDIwmJcPZxSRPJTtAPklsNLMlplZHLgb2DCmzQbgw8H2ncDzPgvWSr9oYRn9g0PsPalhLBHJD1NeymQ63D1pZvcD3wWiwFfcfauZPQRscvcNwJeBx82sCWglFTIAmNk+oByIm9kdwNvdfVs2/xsmctHCcgB2HDvNhQvKQq5GRCTzshogAOMt/+7uD47a7iX1fJHxvrcxo8VNwwW1pcQixs5jHXD5orDLERHJuDk3iR6WeCzC8toSdhzVRLqI5AcFyAxatbBcV2KJSN5QgMygixaWcbith45erYklInOfAmQGXbQwNXm+S70QEckDCpAZtCoIEA1jiUg+UIDMoPrKIsoKY+xUgIhIHlCAzCAzC5Y06Qi7FBGRjFOAzLDhNbFmwc3zIiLTogCZYRfVlXO6N8nR9pxbdV5EZEYpQGbYRWcm0jWMJSJzmwJkhg2vg6UrsURkrlOAzLCKogLqK4vYriVNRGSOU4BkwGWLK3j54KmwyxARySgFSAasbajiYGsPJ05rIl1E5i4FSAasXVoJwJb9bSFXIiKSOQqQDLh4UQXxaISXDmgYS0TmLgVIBiQKolxcX84WBYiIzGEKkAxZ21DFK4fa6U8OhV2KiEhGKEAy5KqlVfQnh9h2VDcUisjcpADJkLUNVQBs2a9hLBGZmxQgGbKwIsGiioTmQURkzlKAZNDapVXqgYjInKUAyaC1DVUcae/lmFbmFZE5SAGSQWuXBvMgGsYSkTlIAZJBa+rKKYxFNIwlInOSAiSD4rEIly+uZOPe1rBLERGZcQqQDLthZQ2vHW7nZGdf2KWIiMwoBUiG3bRqPgA/2dUcciUiIjNLAZJhFy8qp6Y0zo92KkBEZG5RgGRYJGLceGEtP3mjmcEhD7scEZEZowDJgptWzaete4BXDun5ICIydyhAsuDGlTVEDA1jicickvUAMbNbzWynmTWZ2QPjvF9oZk8G7280s8ZR730yOL7TzN6Rzbqno7I4zpUNVfx454mwSxERmTFZDRAziwKPALcBa4B7zGzNmGb3AqfcfQXwMPDZ4HvXAHcDFwO3Al8Mft6scNOFtbxySJfzisjcke0eyNVAk7vvcfd+4Ang9jFtbgceC7afAtabmQXHn3D3PnffCzQFP29W0OW8IjLXxLJ8vnrg4Kj9Q8A1E7Vx96SZtQPVwfEXx3xvfeZKnVmjL+d939rFYZcjIjlqaMhJDjlD7gwOOYPuDA2NbLuT2h4Ktj3VduhMWxhypzxRQEN1cUZrzXaAZJyZ3QfcB9DQ0BByNSMiEePmi+bz3GvH6B0YJFEwa0bfRHKWuzMw6PQPDtGfHGIgeO0btT0wODTqfWdgcCj4Ons7OXp7aIjkoJ/ZHhh0BoNjyaGR95PBL/vh9waHzn5/cNQv/uRg6hd9csjPhMTg0NltZvJS/3ddVscX3r92xn7eeLIdIIeBJaP2FwfHxmtzyMxiQAXQMsnvxd0fBR4FWLduXU7deHHHFfV8c9Mhvr/tOO++fFHY5YjMOHenLzlEd/8g3f1JegcG6e4fpKd/kO6BQXr7B+kZCL76B+lLDtE7MEhvcKxvYIje4FjfqNe+gcEzwdCXDI4lU6GQCdGIEYsYBdEIsagRi0SIRYxYNHVs+P2x7yUKokQiRkHEiI76ikWMSPAajRgRS/2ciBnRCEQjkdSrjbSLRIyojbSPBsciNtJu+PsjNqqNpdrUVRRl5LMZLdsB8ktgpZktI/XL/27g/WPabAA+DPwcuBN43t3dzDYA/2BmnwMWASuBX2St8hlw7fJqFlUkeHrLIQWI5IzhX/odPQN09Cbp6B2gszfJ6d4knX0DwWuSzt4kXf1JOvsG6epLHevuT9LdN0jXqNep/hEdMUgURFNfsQiJgijx4LUwFqGiqIDCssLU8WiEeCxCYSxCYUGEwmA/HosE76W+tyBqFMYiFERTX/Ez2zayHYlQELMz28PhMPzLW84vqwESzGncD3wXiAJfcfetZvYQsMndNwBfBh43syaglVTIELT7JrANSAIfc/fBbNY/XZGIcceV9fz9T/Zw4nQv88sSYZckc4i709Gb5FRXP63d/bR193Oqa4C2ngHau/tp6xmgrXuA9p7UV0fvQCo0epL0D57/L/migiiliRilhTFKCqOUxGPML0tQXJ3aLi6MUhyPUhyPBa9RiuIxiguiFMVTX8XxKIlYajsVGqlf/KnrZGS2MfecGuWZUevWrfNNmzaFXcZZmk6c5m2f+wl/9M7V/PYNy8MuR3LcwOAQJzv7aD498tXS1X/mtaWzj9auflq6+jnV1U9ygj//zaA8UUBlcQEVRamv8uHXRAHlRTHKEwWUJUZeh8OiLFFASTxKLKr7jvOFmW1293XnazfnJtFz3Yr5ZVy2uIKntxxWgOQxd6ejJ8nRjh6OBo89Ptrey/H2Xo519HK8o5fm0320dvcz3t94pYUxqkvjVJfEWTKvmCuWVDKvJM68kjhVxanXyuICqopT+2WJmIZlZMYpQELwvivr+dS3t7H9aAer68rDLkcywN052dnPwVPdHGzt5tCpHg639XA4eD3a1kNX/9kjsBGDmtJCFpQnWFxVxNqlVcwvK6S2rJD5ZQlqywqpKY1TU1qoq/gkJyhAQvDuyxfx6X/ezjMvHVaAzGJDQ87Rjl72nexi78ku9p3sYn9rNwdaujnQ2k3PwNkBMa8kTn1lERfUlnDjyloWVSaoqyiirjJBXUWC2tJCDRPJrKIACUF1aSE3rZrP01sO8/FbLtRfkzmud2CQvSe7eONEJ00nOtl9opPdzZ3sa+mid2Bk8rkwFmFpdTEN80q4fkUNDfOKWDKvmCXziqmvLKKkUP/cZG7R/9Eh+cj1jbz/Sxt5avMhPnjt0rDLEVI9iv2t3ew42sH2Y6fZeayDN46ngmJ4btoMllQVc0FtKiSW15awrLqExpoSFpYnNM8geUUBEpLrLqjmiiWV/N2Pd3P3m5Zo6CLL+pKDvHG8k9cPt/P6kXa2Helgx7HTdAfzEhGDxuoSVi4o5Z2X1bFyQRkr55eyrKZEPUaRgAIkJGbGx966gt/5P5v49qtHeO+VWh8rUwaHnKYTnbx88BSvHmrn1UPt7DjWwcBgqltRloixuq6c31y3hDV15VxUV8bK+WUUxRUUIueiAAnR+ovms2pBGV98YTe3X16v4Y8Z0tbdz0sH2ti8/xRbDqRCo7MvCUBZYYxLF1dw75uXc2l9BZfUl9Mwr1g3somkQQESokjE+OhNF/B7T77MD7Yf5+0XLwy7pFnpaHsPv9jbyi/2tvLLfa3sOt4JpNYzWl1XxnuvrOeKJZVc0VDJsuoSBbXIDFGAhOxdl9Xxv76/ky+80MTbVi/QL7dJONHRy7/tbuHnu1t4cW8L+1u6gdTNdWuXVvHuyxZxVWMVly+u1JVPIhmkf10hi0Uj/OebV/L7T73KtzYf5N+9KXeWoM8VnX1JXtzdwk/faOZfd7fQdCLVwyhPxLhmeTUfuq6Ra5bNY3VdOVEFsEjWKEBywJ1XLeZbmw/xmWd3sH71AmpKC8MuKVRDQ862ox38eFczP97VzJb9p0gOOYmCCFcvq+auqxZz/YoaBYZIyBQgOcDM+Mx7L+XXP/9T/uQ72/j83VeGXVLWdfQO8LM3TvL8jhP8aGfzmWfHX1Jfzu/cuJwbVtZw1dIqCmO6MkokVyhAcsSK+aV89KYL+PwP3+B9axfzlgtrwy4p4w60dPOD7cf54Y7jbNzTSnLIKU/EeMuq+dx0YS03XlhLbVl+98ZEcpkCJId89KYL+PYrR/ijf3qN79x/AxXFBWGXNKPcndcOt/O9rcf5/rbj7Dx+GkiF5703LGP9RQtY21CpmypFZgkFSA5JFET5y7su4+5HX+SjX9/MYx+5moJZ/ss0OTjEL/a18r2tx/nu1mMcbe8lYvCmxnn80TtXc8uaBSytLgm7TBFJgwIkx1y1dB5/9r7L+MS3XuHB//c6n3nvpbPuJre+5CD/1tTCc68f5fvbjnOqe4DCWIQbVtby8VsuZP3qBcwriYddpohMkwIkB9151WL2NHfyxR/tZnlNKb9zY+4/eKqnf5Af7zrBc68f4/ntJzjdl6SsMMb61fN5x8ULecuqWorj+t9NZC7Rv+gc9Ym3r2JfSxd/+ux2Bt353RuX51xPpKN3gBd2nOC5147xo10n6B0Yoqq4gF+/tI5bL13I9RfUEI/N7iE4EZmYAiRHRSLG537zCsxe4c+f28He5i7+5I5LQv+F3Hy6jx9sP86/vH6Mf9t9koFBZ35ZIXddtYTbLlnI1cvmaRJcJE8oQHJYoiDK39x9JctrSvib55s40NrNX951GYurirNWg3tqJdsfbD/B97cd46WDbbhDw7xifuvXGnnHxQtZ21ClJVhE8pC5e9g1ZMy6det806ZNYZcxI57ecoj/9sxruKcu9/3dGy/I2HLjnX1JNu5p4Uc7m3lh5wkOneoBUjf13bJ6IbesWcDqurKcG1ITkZlhZpvdfd152ylAZo/DbT382bPb+c6rR1lUkeAD1y7lvVfWs6iyaFo/t6N3ILX8+b5Wfr6nhZcOtJEccooKoly/opq3XjSft66aP+3ziMjsoABh7gXIsBf3tPC57+/iF3tbMYPrlldz3fJqLllcwSWLKqgpjY/bO+jsS3L4VA+H27p543gn2452sP1oB2+c6MQ99RS+S+sruH5FDW9eUcPapVV6+p5IHlKAMHcDZNiBlm6efukQ//zqUZqaUyEAqedglCdilBcVMOROT/8QPf1JuoLHtQ6rq0iwuq6cS+sreFPjPK5oqKRUy5+L5D0FCHM/QEY73TvA1iMdbDvSQUtXH+09A3T0JIlFjEQ8SiIWpbaskPqqIuori1heU0KVbuYTkXFMNkD05+YcUZYo4Nrl1Vy7vDrsUkQkT+iCfRERSYsCRERE0qIAERGRtChAREQkLQoQERFJiwJERETSogAREZG0KEBERCQtc/pOdDNrBvaHXcc01QAnwy4ih+jzOJs+jxH6LM42nc9jqbvXnq/RnA6QucDMNk1mSYF8oc/jbPo8RuizOFs2Pg8NYYmISFoUICIikhYFSO57NOwCcow+j7Pp8xihz+JsGf88NAciIiJpUQ9ERETSogDJcWYWNbOXzOw7YdcSNjPbZ2avmdnLZpYfTwqbgJlVmtlTZrbDzLab2XVh1xQWM1sV/D8x/NVhZr8Xdl1hMbP/YmZbzex1M/uGmSUydi4NYeU2M/s4sA4od/d3hV1PmMxsH7DO3fP+Wn8zewz4qbt/ycziQLG7t4VdV9jMLAocBq5x99l+D9iUmVk98DNgjbv3mNk3gWfd/WuZOJ96IDnMzBYD7wS+FHYtkjvMrAK4EfgygLv3KzzOWA/szsfwGCUGFJlZDCgGjmTqRAqQ3Pa/gT8AhsIuJEc48D0z22xm94VdTIiWAc3AV4PhzS+ZWUnYReWIu4FvhF1EWNz9MPBXwAHgKNDu7t/L1PkUIDnKzN4FnHD3zWHXkkPe7O5rgduAj5nZjWEXFJIYsBb4W3e/EugCHgi3pPAFQ3nvAb4Vdi1hMbMq4HZSf2QsAkrM7IOZOp8CJHddD7wnGPd/ArjZzP5vuCWFK/jrCnc/ATwDXB1uRaE5BBxy943B/lOkAiXf3QZscffjYRcSorcBe9292d0HgKeBX8vUyRQgOcrdP+nui929kVS3/Hl3z9hfErnOzErMrGx4G3g78Hq4VYXD3Y8BB81sVXBoPbAtxJJyxT3k8fBV4ABwrZkVm5mR+n9je6ZOFsvUDxaZYQuAZ1L/JogB/+Du/xJuSaH6T8DXg2GbPcC/D7meUAV/VNwC/G7YtYTJ3Tea2VPAFiAJvEQG70jXZbwiIpIWDWGJiEhaFCAiIpIWBYiIiKRFASIiImlRgIiISFoUICIikhYFiIiIpEXa7WRrAAABYUlEQVQBIhIyM/vt4DkneX0zoMw+ChCR8P0GcDNwV9iFiEyFAkQkBGb2d2Z2fbC7ETgRvIrMGgoQkSwKnpgHcC3wYrBdCvwUqAilKJE0aTFFkQwzs28BrcDlwHfM7B+BXe4+aGYR4H2khrBeMLPfd/fBEMsVmTT1QEQy71LguLtf6+6fJvXciuGVhG8GXnH3fcArwb7IrKAAEckgM0sA84CHRh1+ByMB8gFGnmHxjWBfZFbQEJZIZl0MbHT3JICZFQOV7n7EzIpIPX50vZn9Bak/6MrMrMjde8IrWWRy1AMRyaxLgVdH7b8VeCHYfg/wnLs3uHujuzcA3wbeneUaRdKiABHJrLEBMnr+4wOknu0+2jNA3j66WGYXPZFQJIvMbAtwjbsPhF2LyHQpQEREJC0awhIRkbQoQEREJC0KEBERSYsCRERE0qIAERGRtChAREQkLQoQERFJiwJERETS8v8B2/RoIeYpP9sAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "def lj_force(r, epsilon, sigma):\n", " \"\"\"\n", " Implementation of the Lennard-Jones potential \n", " to calculate the force of the interaction.\n", " \n", " Parameters\n", " ----------\n", " r: float\n", " Distance between two particles (Å)\n", " epsilon: float \n", " Potential energy at the equilibrium bond \n", " length (eV)\n", " sigma: float \n", " Distance at which the potential energy is \n", " zero (Å)\n", " \n", " Returns\n", " -------\n", " float\n", " Force of the van der Waals interaction (eV/Å)\n", " \"\"\"\n", " return 48 * epsilon * np.power(\n", " sigma, 12) / np.power(\n", " r, 13) - 24 * epsilon * np.power(\n", " sigma, 6) / np.power(r, 7)\n", " \n", "r = np.linspace(3.5, 8, 100)\n", "plt.plot(r, lj_force(r, 0.0103, 3.4))\n", "plt.xlabel(r'$r$/Å')\n", "plt.ylabel(r'$f$/eVÅ$^{-1}$')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You may have noticed that in Newton's second law of motion, the force is a vector quantity, whereas the first negative derivative of the energy is a scalar.\n", "Therefore, it is important that we determine the force in each direction for our simulation.\n", "This is achieved by multiplication by the unit vector in that direction,\n", "\n", "$$ \\mathbf{f}_x = f\\hat{\\mathbf{r}}_x\\text{, where }\\hat{\\mathbf{r}}_x = \\dfrac{r_x}{|\\mathbf{r}|}. $$\n", "\n", "In the above equation, $r_x$ is the distance between the two particles in the $x$-dimension and $|\\mathbf{r}|$ is the magnitude of the distance vector. \n", "For simplicity, we will initially **only** consider particles interacting in a one-dimensional space. \n", "\n", "The Python code below shows how to determine the acceleration on each atom of argon due to each other atom of argon. \n", "It is possible to **increase the efficiency** of this algorithm by applying Newton's third law, e.g. the force on atom $i$ will be equal and opposite to the force on atom $j$. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Acceleration on particle 0 = 1.453e-04 eV/Åamu\n", "Acceleration on particle 1 = -4.519e-05 eV/Åamu\n", "Acceleration on particle 2 = -1.002e-04 eV/Åamu\n" ] } ], "source": [ "mass_of_argon = 39.948 # amu\n", "\n", "def get_accelerations(positions):\n", " \"\"\"\n", " Calculate the acceleration on each \n", " particle as a result of each other \n", " particle. \n", " N.B. We use the Python convention of\n", " numbering from 0.\n", " \n", " Parameters\n", " ----------\n", " positions: ndarray of floats\n", " The positions, in a single dimension, \n", " for all of the particles (Å)\n", " \n", " Returns\n", " -------\n", " ndarray of floats\n", " The acceleration on each particle (eV/Åamu)\n", " \"\"\"\n", " accel_x = np.zeros((positions.size, positions.size))\n", " for i in range(0, positions.size - 1):\n", " for j in range(i + 1, positions.size):\n", " r_x = positions[j] - positions[i]\n", " rmag = np.sqrt(r_x * r_x)\n", " force_scalar = lj_force(rmag, 0.0103, 3.4)\n", " force_x = force_scalar * r_x / rmag\n", " accel_x[i, j] = force_x / mass_of_argon #eV Å-1 amu-1\n", " # appling Newton's third law\n", " accel_x[j, i] = - force_x / mass_of_argon\n", " return np.sum(accel_x, axis=0)\n", "\n", "accel = get_accelerations(np.array([1, 5, 10]))\n", "print('Acceleration on particle 0 = {:.3e} eV/Åamu'.format(\n", " accel[0]))\n", "print('Acceleration on particle 1 = {:.3e} eV/Åamu'.format(\n", " accel[1]))\n", "print('Acceleration on particle 2 = {:.3e} eV/Åamu'.format(\n", " accel[2]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Integration\n", "\n", "Now that we have seen how to obtain the acceleration on our particles, we can apply the **Newtonian equations of motion** to probe the particles trajectory, \n", "\n", "$$ \\mathbf{x}_i(t + \\Delta t) = \\mathbf{x}_i(t) + \\mathbf{v}_i(t)\\Delta t + \\dfrac{1}{2} \\mathbf{a}_i(t)\\Delta t^2, $$\n", "\n", "$$ \\mathbf{v}_i(t + \\Delta t) = \\mathbf{v}_i(t) + \\dfrac{1}{2}\\big[\\mathbf{a}_i(t) + \\mathbf{a}_i(t+\\Delta t)\\big]\\Delta t, $$ \n", "\n", "where, $\\Delta t$ is the timestep (how far in time is incremented), $\\mathbf{x}_i$ is the particle position, $\\mathbf{v}_i$ is the velocity, and $\\mathbf{a}_i$ the acceleration. \n", "This pair of equations is known as the Velocity-Verlet algorithm, which can be written as:\n", "\n", "1. Calculate the force (and therefore acceleration) on the particle\n", "2. Find the position of the particle after some timestep\n", "3. Calculate the new forces and accelerations\n", "4. Determine a new velocity for the particle, based on the average acceleration at the current and new positions\n", "5. Overwrite the old acceleration values with the new ones, $\\mathbf{a}_i(t) = \\mathbf{a}_i(t+\\Delta t)$\n", "6. Repeat\n", "\n", "After the initial relaxation of the particles to equilibrium, this process can be continued for as long as is required to get sufficient statistics for the quantity you are intereseting in. \n", "\n", "The Python code below is a set of two function for the above equations, these will be applied later." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def update_pos(x, v, a, dt):\n", " \"\"\"\n", " Update the particle positions.\n", " \n", " Parameters\n", " ----------\n", " x: ndarray of floats\n", " The positions of the particles in a \n", " single dimension (Å)\n", " v: ndarray of floats\n", " The velocities of the particles in a \n", " single dimension (eVs/Åamu)\n", " a: ndarray of floats\n", " The accelerations of the particles in a \n", " single dimension (eV/Åamu)\n", " dt: float\n", " The timestep length (s)\n", " \n", " Returns\n", " -------\n", " ndarray of floats:\n", " New positions of the particles in a \n", " single dimension (Å)\n", " \"\"\"\n", " return x + v * dt + 0.5 * a * dt * dt\n", "\n", "def update_velo(v, a, a1, dt):\n", " \"\"\"\n", " Update the particle velocities.\n", " \n", " Parameters\n", " ----------\n", " v: ndarray of floats\n", " The velocities of the particles in a \n", " single dimension (eVs/Åamu)\n", " a: ndarray of floats\n", " The accelerations of the particles in a \n", " single dimension at the previous \n", " timestep (eV/Åamu)\n", " a1: ndarray of floats\n", " The accelerations of the particles in a \n", " single dimension at the current \n", " timestep (eV/Åamu)\n", " dt: float\n", " The timestep length (s)\n", " \n", " Returns\n", " -------\n", " ndarray of floats:\n", " New velocities of the particles in a \n", " single dimension (eVs/Åamu)\n", " \"\"\"\n", " return v + 0.5 * (a + a1) * dt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above process is called the intergration step, and the Velocity-Verlet algorithm is the **integrator**. \n", "This function is highly non-linear for more than two particles. \n", "The result is that the integration step will only be valid for very small values of $\\Delta t$, e.g. if a large timestep is used the acceleration calculated will not be accurate as the forces on the atom will change too significantly during it. \n", "The value for the timestep is usually on the order of 10-15 s (femtoseconds). \n", "So in order to measure a nanosecond of \"real-time\" molecular dynamics, 1 000 000 (one million) iterations of the above algorithm must be performed. \n", "This can be very slow for large, realistic systems.\n", "\n", "## Initialisation\n", "\n", "There are only **two components** left that we need to run a molecular dynamics simulation, and both are associated with the original configuration of the system; the original particle positions, and the original particle velocities. \n", "\n", "The particle positions are usually either taken from some library of structures (e.g. the [protein data bank](http://www.rcsb.org) if you are simulating proteins) or based on some knowledge of the system (e.g. CaF2 is well known to have a face-centred cubic structure). \n", "For complex, multicomponent systems, software such as Packmol [[1](#References)] may be used to build up the structure from its constituent parts.\n", "The importance of this initial structure **cannot be overstated**.\n", "For example if the initial structure is unrepresentative of the equilibrium structure, it may take a very long time before the equilibrium structure is obtained, possibly much longer than can be reasonably simulated. \n", "\n", "The particle velocities are more general, as the total kinetic energy, $E_K$ of the system (and therefore the particle velocities) are dependent on the temperature of the simulation, $T$. \n", "\n", "$$ E_K = \\sum_{i=1}^N \\dfrac{m_i|v_i|^2}{2} = \\dfrac{3}{2}Nk_BT, $$\n", "\n", "where, $m_i$ is the mass of particle $i$, $N$ is the number of particles, and $k_B$ is the Boltzmann constant. \n", "A common method to initialise the velocities is implemented in the Python function below. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "from scipy.constants import Boltzmann\n", "\n", "def init_velocity(T, number_of_particles):\n", " \"\"\"\n", " Initialise the velocities for a series \n", " of particles.\n", " \n", " Parameters\n", " ----------\n", " T: float\n", " Temperature of the system at \n", " initialisation (K)\n", " number_of_particles: int\n", " Number of particles in the \n", " system\n", " \n", " Returns\n", " -------\n", " ndarray of floats\n", " Initial velocities for a series of \n", " particles (eVs/Åamu)\n", " \"\"\"\n", " R = np.random.rand(number_of_particles) - 0.5\n", " return R * np.sqrt(Boltzmann * T / (mass_of_argon * 1.602e-19))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "\n", "1. Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. *J. Comput. Chem.* 2009, **30** (13), 2157–2164. [10.1002/jcc.21224](https://doi.org/10.1002/jcc.21224).\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Python 3", "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.7.1" } }, "nbformat": 4, "nbformat_minor": 2 }