{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Molecular dynamics simulations of bulk water\n", "\n", "In this example, we show how to perform molecular dynamics of bulk water using the popular interatomic TIP3P potential\n", "([W. L. Jorgensen et. al.](https://doi.org/10.1063/1.445869)) and [LAMMPS](http://lammps.sandia.gov/)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "%matplotlib inline \n", "import matplotlib.pylab as plt\n", "from pyiron.project import Project\n", "import ase.units as units" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "pr = Project(\"tip3p_water\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating the initial structure\n", "\n", "We will setup a cubic simulation box consisting of 27 water molecules density density is 1 g/cm\$^3\$. The target density is achieved by determining the required size of the simulation cell and repeating it in all three spatial dimensions" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "2635f5900e0d4ab29f46adec6f56e2a7", "version_major": 2, "version_minor": 0 }, "text/html": [ "

Failed to display Jupyter Widget of type `NGLWidget`.

\n", "

\n", " If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n", " that the widgets JavaScript is still loading. If this message persists, it\n", " likely means that the widgets JavaScript library is either not installed or\n", " not enabled. See the Jupyter\n", " Widgets Documentation for setup instructions.\n", "

\n", "

\n", " If you're reading this message in another frontend (for example, a static\n", " rendering on GitHub or NBViewer),\n", " it may mean that your frontend doesn't currently support widgets.\n", "

\n" ], "text/plain": [ "NGLWidget()" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "density = 1.0e-24 # g/A^3\n", "n_mols = 27\n", "mol_mass_water = 18.015 # g/mol\n", "\n", "# Determining the supercell size size\n", "mass = mol_mass_water * n_mols / units.mol # g\n", "vol_h2o = mass / density # in A^3\n", "a = vol_h2o ** (1./3.) # A\n", "\n", "# Constructing the unitcell\n", "n = int(round(n_mols ** (1. / 3.)))\n", "\n", "dx = 0.7\n", "r_O = [0, 0, 0]\n", "r_H1 = [dx, dx, 0]\n", "r_H2 = [-dx, dx, 0]\n", "unit_cell = (a / n) * np.eye(3)\n", "water = pr.create_atoms(elements=['H', 'H', 'O'], \n", " positions=[r_H1, r_H2, r_O], \n", " cell=unit_cell)\n", "water.set_repeat([n, n, n])\n", "water.plot3d()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'H54O27'" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "water.get_chemical_formula()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Equilibrate water structure\n", "\n", "The initial water structure is obviously a poor starting point and requires equilibration (Due to the highly artificial structure a MD simulation with a standard time step of 1fs shows poor convergence). Molecular dynamics using a time step that is two orders of magnitude smaller allows us to generate an equilibrated water structure. We use the NVT ensemble for this calculation:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['H2O_tip3p', 'H2O_tip3p_slab']" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "job_name = \"water_slow\"\n", "ham = pr.create_job(\"Lammps\", job_name)\n", "ham.structure = water\n", "# Listing available potentials for this structure\n", "ham.list_potentials()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will use the `H2O_tip3p` potential. The `H2O_tip3p_slab` should be used if you want to switch off the periodic boundary conditions along the z-axis " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "ham.potential = 'H2O_tip3p'\n", "ham.calc_md(temperature=300, \n", " n_ionic_steps=1e4, \n", " time_step=0.01)\n", "ham.run()" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "scrolled": true }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "383b3d53a6324a9eace3d6bb81c16f58", "version_major": 2, "version_minor": 0 }, "text/html": [ "

Failed to display Jupyter Widget of type `NGLWidget`.

\n", "

\n", " If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n", " that the widgets JavaScript is still loading. If this message persists, it\n", " likely means that the widgets JavaScript library is either not installed or\n", " not enabled. See the Jupyter\n", " Widgets Documentation for setup instructions.\n", "

\n", "

\n", " If you're reading this message in another frontend (for example, a static\n", " rendering on GitHub or NBViewer),\n", " it may mean that your frontend doesn't currently support widgets.\n", "

\n" ], "text/plain": [ "NGLWidget(count=101)" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "view = ham.animate_structure()\n", "view" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Full equilibration\n", "\n", "At the end of this simulation, we have obtained a structure that approximately resembles water. Now we increase the time step to get a reasonably equilibrated structure " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Get the final structure from the previous simulation\n", "struct = ham.get_structure(iteration_step=-1)\n", "job_name = \"water_fast\"\n", "ham_eq = pr.create_job(\"Lammps\", job_name)\n", "ham_eq.structure = struct\n", "ham_eq.potential = 'H2O_tip3p'\n", "ham_eq.calc_md(temperature=300, \n", " n_ionic_steps=1e4,\n", " n_print=10, \n", " dt=1)\n", "ham_eq.run()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "e9996981cf4047d4bca6bd6a62db4fb0", "version_major": 2, "version_minor": 0 }, "text/html": [ "

Failed to display Jupyter Widget of type `NGLWidget`.

\n", "

\n", " If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n", " that the widgets JavaScript is still loading. If this message persists, it\n", " likely means that the widgets JavaScript library is either not installed or\n", " not enabled. See the Jupyter\n", " Widgets Documentation for setup instructions.\n", "

\n", "

\n", " If you're reading this message in another frontend (for example, a static\n", " rendering on GitHub or NBViewer),\n", " it may mean that your frontend doesn't currently support widgets.\n", "