From 3faca27d4c1704dfe68211d498001871977bd608 Mon Sep 17 00:00:00 2001 From: Levi Petix <81758680+clpetix@users.noreply.github.com> Date: Mon, 8 Jul 2024 06:39:34 -0500 Subject: [PATCH] Add Multiparticle Collision Dynamics tutorial (#122) * Add MPCD tutorial. * Update tutorial to not rely on external gsd file. * Update authors list. * Delete .DS_Store * Update Diffusion in Solution tutorial. * Delete 09-Multiparticle-Collision-Dynamics/.DS_Store * Add semicolons to supress output. * Include MPCD in outline. * Test with `trunk-minor`. * Ruff fixes and update matplotlib imports and usage. * Update pressure driven flow example. * Update diffusion in solution example. * Fix mpcd_particle_sorter in Diffusion in Solution example. * Update to use Frame rather than Snapshot. --------- Co-authored-by: Joshua A. Anderson --- .../00-index.ipynb | 65 + .../01-Pressure-Driven-Flow.ipynb | 1694 +++++++++++++++++ .../02-Diffusion-in-Solution.ipynb | 1364 +++++++++++++ AUTHORS.md | 3 + README.md | 1 + 5 files changed, 3127 insertions(+) create mode 100644 09-Multiparticle-Collision-Dynamics/00-index.ipynb create mode 100644 09-Multiparticle-Collision-Dynamics/01-Pressure-Driven-Flow.ipynb create mode 100644 09-Multiparticle-Collision-Dynamics/02-Diffusion-in-Solution.ipynb diff --git a/09-Multiparticle-Collision-Dynamics/00-index.ipynb b/09-Multiparticle-Collision-Dynamics/00-index.ipynb new file mode 100644 index 0000000..e00c9ba --- /dev/null +++ b/09-Multiparticle-Collision-Dynamics/00-index.ipynb @@ -0,0 +1,65 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Introducing Multiparticle Collison Dynamics\n", + "\n", + "This tutorial introduces how to use the multiparticle collision dynamics (MPCD)\n", + "method in HOOMD-blue with examples of pressure-driven flow and diffusion of a\n", + "solution of nearly-hard particles.\n", + "\n", + "**Prerequisites:**\n", + "\n", + "This tutorial assumes you are familiar with molecular dynamics (MD) (see\n", + "[Introducing Molecular Dynamics](../01-Introducing-Molecular-Dynamics/00-index.ipynb))\n", + "and have background knowledge on MPCD. To learn more about MPCD, consider\n", + "reading one of the following review articles:\n", + "\n", + "* [Modeling hydrodynamic interactions in soft materials with multiparticle collision dynamics](https://doi.org/10.1016/j.coche.2019.02.007)\n", + "* [Multi-Particle Collision Dynamics: A Particle-Based Mesoscale Simulation Approach to the Hydrodynamics of Complex Fluids](https://doi.org/10.1007/978-3-540-87706-6_1)\n", + "* [Multiparticle Collision Dynamics: Simulation of Complex Systems on Mesoscales](https://doi.org/10.1002/9780470371572.ch2)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "nbsphinx-toctree": { + "maxdepth": 1 + } + }, + "source": [ + "## Outline\n", + "\n", + "1. [Pressure-Driven Flow Between Parallel Plates](01-Pressure-Driven-Flow.ipynb)\n", + "2. [Diffusion of a Solution of Nearly-Hard Spheres](02-Diffusion.ipynb)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This tutorial is written with [jupyter](https://jupyter.org/). You can download the source from the [hoomd-examples](https://github.com/glotzerlab/hoomd-examples) repository." + ] + } + ], + "metadata": { + "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.12.2" + }, + "record_timing": false + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/09-Multiparticle-Collision-Dynamics/01-Pressure-Driven-Flow.ipynb b/09-Multiparticle-Collision-Dynamics/01-Pressure-Driven-Flow.ipynb new file mode 100644 index 0000000..8e68ecd --- /dev/null +++ b/09-Multiparticle-Collision-Dynamics/01-Pressure-Driven-Flow.ipynb @@ -0,0 +1,1694 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0", + "metadata": {}, + "source": [ + "# Pressure-Driven Flow Between Parallel Plates" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, + "source": [ + "## Overview" + ] + }, + { + "cell_type": "markdown", + "id": "2", + "metadata": {}, + "source": [ + "### Questions" + ] + }, + { + "cell_type": "markdown", + "id": "3", + "metadata": {}, + "source": [ + "* What is multiparticle collision dynamics (MPCD)?\n", + "* How do I set up a pressure-driven flow simulation using MPCD?" + ] + }, + { + "cell_type": "markdown", + "id": "4", + "metadata": {}, + "source": [ + "### Objectives" + ] + }, + { + "cell_type": "markdown", + "id": "5", + "metadata": {}, + "source": [ + "* Demonstrate how to initialize MPCD particles using a **Snapshot**.\n", + "* Explain how to configure the MPCD **Integrator** with a **StreamingMethod**\n", + " and **CollisionMethod**.\n", + "* Show how to include a confining **Geometry** in both the streaming and\n", + " collision steps.\n", + "* Show how to apply a **BodyForce** to MPCD particles.\n", + "* Demonstrate how to run a simulation and gather data." + ] + }, + { + "cell_type": "markdown", + "id": "6", + "metadata": {}, + "source": [ + "## Boilerplate code" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "7", + "metadata": {}, + "outputs": [], + "source": [ + "import hoomd\n", + "import matplotlib\n", + "import numpy\n", + "\n", + "%matplotlib inline\n", + "matplotlib.style.use('ggplot')\n", + "import matplotlib_inline\n", + "\n", + "matplotlib_inline.backend_inline.set_matplotlib_formats('svg')" + ] + }, + { + "cell_type": "markdown", + "id": "8", + "metadata": {}, + "source": [ + "## Initialization" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, + "source": [ + "MPCD particles are currently initialized using a `Snapshot`. We will place the\n", + "particles between two parallel plates. The parallel plates have a surface normal\n", + "in the *y* direction, and we have chosen the separation distance between them as\n", + "$2 H = 20$. The plates are infinite in the *x* and *z* directions, so\n", + "periodic boundaries will be used in these directions. We choose the extent in\n", + "the periodic directions to be $L = 25$.\n", + "\n", + "The simulation box needs padded in the *y* direction to account for the\n", + "collision cells wrapping through the *y* periodic boundary. A padding of 4\n", + "cells (2 each direction) should be sufficient." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "10", + "metadata": {}, + "outputs": [], + "source": [ + "H = 10\n", + "L = 25\n", + "padding = 4\n", + "snapshot = hoomd.Snapshot()\n", + "snapshot.configuration.box = [L, 2*H + padding, L, 0, 0, 0]" + ] + }, + { + "cell_type": "markdown", + "id": "11", + "metadata": {}, + "source": [ + "The MPCD particle positions will be randomly drawn between the parallel plates\n", + "at a target number density $\\rho = 5$. All particles will use the default\n", + "typeid (0), which we will call \"A\", and default mass (1). The number of\n", + "particles *N* can be computed using the density, $N = \\rho V$ where $V = 2 H\n", + "L^2$ is the volume between the parallel plates." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "12", + "metadata": {}, + "outputs": [], + "source": [ + "density = 5\n", + "rng = numpy.random.default_rng(seed=42)\n", + "\n", + "snapshot.mpcd.types = [\"A\"]\n", + "snapshot.mpcd.N = int(density * (2*H) * L**2)\n", + "snapshot.mpcd.position[:] = [0.5*L, H, 0.5*L] * rng.uniform(\n", + " low=-1, high=1, size=(snapshot.mpcd.N, 3))" + ] + }, + { + "cell_type": "markdown", + "id": "13", + "metadata": {}, + "source": [ + "The particle velocities also need to be initialized for the thermostat to work\n", + "properly with the stochastic rotation dynamics (SRD) collision rule (see below).\n", + "We randomly draw them to be consistent with the Maxwell-Boltzmann distribution\n", + "for $k_{\\rm B} T = 1$ and unit mass, and we force them to have zero mean." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "14", + "metadata": {}, + "outputs": [], + "source": [ + "kT = 1\n", + "velocity = rng.normal(0.0, numpy.sqrt(kT), (snapshot.mpcd.N, 3))\n", + "velocity -= numpy.mean(velocity, axis=0)\n", + "snapshot.mpcd.velocity[:] = velocity" + ] + }, + { + "cell_type": "markdown", + "id": "15", + "metadata": {}, + "source": [ + "The `Snapshot` is now fully initialized, so we can use it to create a `Simulation`.\n", + "Note that this works even though there are no regular HOOMD particles in the\n", + "snapshot." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "16", + "metadata": {}, + "outputs": [], + "source": [ + "simulation = hoomd.Simulation(device=hoomd.device.CPU(), seed=1)\n", + "simulation.create_state_from_snapshot(snapshot)" + ] + }, + { + "cell_type": "markdown", + "id": "17", + "metadata": {}, + "source": [ + "## Configuring the MPCD integrator" + ] + }, + { + "cell_type": "markdown", + "id": "18", + "metadata": {}, + "source": [ + "The motion of the MPCD particles is governed by alternating streaming and\n", + "collision steps. During the streaming step, the particles move according to\n", + "Newton's equations of motion. During the collision step, the particles are\n", + "assigned to a \"collision\" cell and stochastically exchange momentum with each\n", + "other in a cell. The time between collisions is called the collision time. These\n", + "equations of motion are implemented by the MPCD `Integrator`. The amount of time\n", + "covered by a streaming step or the time between collisions is a multiple of the\n", + "timestep of this integrator. Refer to the `Integrator` documentation for more\n", + "information about when and the order specific steps of the algorithm trigger.\n", + "\n", + "Here we will use a timestep of $\\Delta t = 0.1$, which will correspond to our\n", + "collision time (see below)." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "19", + "metadata": {}, + "outputs": [], + "source": [ + "integrator = hoomd.mpcd.Integrator(dt=0.1)" + ] + }, + { + "cell_type": "markdown", + "id": "20", + "metadata": {}, + "source": [ + "We will use the SRD collision rule with collision time $\\Delta t = 0.1$,\n", + "collision angle $\\alpha = 130^\\circ$, and a thermostat to maintain constant\n", + "temperature. We attach an SRD collision method to the `Integrator` and have a\n", + "collision occur every timestep." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "21", + "metadata": {}, + "outputs": [], + "source": [ + "integrator.collision_method = hoomd.mpcd.collide.StochasticRotationDynamics(\n", + " period=1, angle=130, kT=kT)" + ] + }, + { + "cell_type": "markdown", + "id": "22", + "metadata": {}, + "source": [ + "Next, we setup the streaming method. The streaming method will use bounce back\n", + "integration between parallel plates, which is defined by a geometry." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "23", + "metadata": {}, + "outputs": [], + "source": [ + "plates = hoomd.mpcd.geometry.ParallelPlates(separation=2*H, speed=0.0, no_slip=True)" + ] + }, + { + "cell_type": "markdown", + "id": "24", + "metadata": {}, + "source": [ + "The streaming method will also apply a constant force $f_x$ in the $x$ direction\n", + "to create a flow. For pressure-driven flow between parallel plates, the\n", + "expected velocity field has the parabolic form\n", + "$$\n", + "v_x(y) = \\frac{\\rho f_x H^2}{2 \\mu}\\left[1 - \\left(\\frac{y}{H}\\right)^2 \\right]\n", + "$$\n", + "where $\\mu = 3.96$ is the dynamic viscosity of the MPCD solvent for the chosen\n", + "parameters. We will use a force $f_x = 0.004$ to give a theoretically predicted maximum\n", + "velocity of 0.25." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "25", + "metadata": {}, + "outputs": [], + "source": [ + "fx = 0.004\n", + "body_force = hoomd.mpcd.force.ConstantForce(force=(fx, 0, 0))" + ] + }, + { + "cell_type": "markdown", + "id": "26", + "metadata": {}, + "source": [ + "We combine all these together to make the streaming method. We stream and\n", + "collide with the same period because a constant force is easily integrated using\n", + "the velocity Verlet method, and better performance is usually achieved when\n", + "fewer streaming steps are needed." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "27", + "metadata": {}, + "outputs": [], + "source": [ + "integrator.streaming_method = hoomd.mpcd.stream.BounceBack(\n", + " period=1,\n", + " geometry=plates,\n", + " mpcd_particle_force=body_force\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "28", + "metadata": {}, + "source": [ + "A virtual particle filler is also needed to ensure the collision method functions\n", + "correctly with the bounce-back geometry. We can attach one to the integrator\n", + "using the same `plates` geometry we defined above." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "29", + "metadata": {}, + "outputs": [], + "source": [ + "filler = hoomd.mpcd.fill.GeometryFiller(\n", + " type=\"A\", density=density, kT=kT, geometry=plates)\n", + "integrator.virtual_particle_fillers.append(filler)" + ] + }, + { + "cell_type": "markdown", + "id": "30", + "metadata": {}, + "source": [ + "The performance of MPCD simulations benefits from periodically sorting the\n", + "particles into cell order. It is recommended to attach a tuner to the integrator\n", + "to do this. The optimal number of collision steps between sorting varies based\n", + "on the simulation and your computing hardware, but typically tens of collisions\n", + "works well." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "31", + "metadata": {}, + "outputs": [], + "source": [ + "integrator.mpcd_particle_sorter = hoomd.mpcd.tune.ParticleSorter(trigger=20)" + ] + }, + { + "cell_type": "markdown", + "id": "32", + "metadata": {}, + "source": [ + "Last, don't forget to attach the integrator to the simulation operations!" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "33", + "metadata": {}, + "outputs": [], + "source": [ + "simulation.operations.integrator = integrator" + ] + }, + { + "cell_type": "markdown", + "id": "34", + "metadata": {}, + "source": [ + "## Run simulation" + ] + }, + { + "cell_type": "markdown", + "id": "35", + "metadata": {}, + "source": [ + "Now that the simulation is configured, we will first run a warmup period to\n", + "allow the system to come to steady state." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "36", + "metadata": {}, + "outputs": [], + "source": [ + "simulation.run(10000)" + ] + }, + { + "cell_type": "markdown", + "id": "37", + "metadata": {}, + "source": [ + "After reaching steady state, we can measure $v_x(y)$ using histogramming. We use bin spacing 0.1." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "38", + "metadata": {}, + "outputs": [], + "source": [ + "bin_size = 0.2\n", + "num_bins = int(2 * H / bin_size)\n", + "bin_edges = numpy.linspace(-H, H, num_bins + 1)\n", + "bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])\n", + "\n", + "num_samples = 0\n", + "counts = numpy.zeros(num_bins, dtype=int)\n", + "total_velocity = numpy.zeros(num_bins, dtype=float)" + ] + }, + { + "cell_type": "markdown", + "id": "39", + "metadata": {}, + "source": [ + "Then, we run the simulation long enough to get enough independent samples with\n", + "good statistics." + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "40", + "metadata": {}, + "outputs": [], + "source": [ + "for _ in range(500):\n", + " simulation.run(100)\n", + "\n", + " # take a snapshot of the system\n", + " snap = simulation.state.get_snapshot()\n", + " y = snap.mpcd.position[:, 1]\n", + " v_x = snap.mpcd.velocity[:,0]\n", + "\n", + " # bin x-velocity using the y-position, and accumulate histograms\n", + " y_bins = numpy.digitize(y, bin_edges) - 1\n", + " counts += numpy.bincount(y_bins, minlength=num_bins)\n", + " total_velocity += numpy.bincount(y_bins, weights=v_x, minlength=num_bins)\n", + " num_samples += 1\n", + "\n", + "average_velocity = numpy.zeros(num_bins, dtype=float)\n", + "numpy.divide(total_velocity, counts, out=average_velocity, where=(counts > 0));" + ] + }, + { + "cell_type": "markdown", + "id": "41", + "metadata": {}, + "source": [ + "Now, we plot our results for the average velocity profile and compare it to\n", + "the theoretical expectation. We see that the agreement is generally good, but\n", + "there is a small amount of slip near the walls. This is typical of MPCD\n", + "simulations." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "42", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-06-18T11:12:07.900188\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.9.0, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "fig = matplotlib.figure.Figure(figsize=(5, 3.09))\n", + "ax = fig.add_subplot()\n", + "\n", + "ax.plot(bin_centers, average_velocity, label=\"simulated\")\n", + "\n", + "viscosity = 3.96\n", + "vx_theory = ((density * fx * H**2) / (2*viscosity)) * (1 - (bin_centers/H)**2)\n", + "ax.plot(bin_centers, vx_theory, ls=\"--\", label=\"expected\")\n", + "\n", + "ax.set_xlabel(r\"$y$\")\n", + "ax.set_xlim([-H, H])\n", + "ax.set_ylabel(r\"$v_x$\")\n", + "ax.set_ylim([0.0, 0.3])\n", + "ax.legend()\n", + "fig" + ] + } + ], + "metadata": { + "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.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/09-Multiparticle-Collision-Dynamics/02-Diffusion-in-Solution.ipynb b/09-Multiparticle-Collision-Dynamics/02-Diffusion-in-Solution.ipynb new file mode 100644 index 0000000..197afea --- /dev/null +++ b/09-Multiparticle-Collision-Dynamics/02-Diffusion-in-Solution.ipynb @@ -0,0 +1,1364 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Diffusion of a Solution of Nearly-Hard Spheres" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Overview" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Questions" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* How do I set up a simulation of of MD particles (solutes) in MPCD particles\n", + " (solvent)?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Objectives" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* Demonstrate how to add MPCD particles to an existing MD particle **Frame**\n", + " or GSD file.\n", + "* Show how to couple MD particles to MPCD particles through the **CollisionMethod**." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Boilerplate Code" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import itertools\n", + "\n", + "import freud\n", + "import gsd.hoomd\n", + "import hoomd\n", + "import matplotlib\n", + "import numpy\n", + "\n", + "%matplotlib inline\n", + "matplotlib.style.use('ggplot')\n", + "import matplotlib_inline\n", + "\n", + "matplotlib_inline.backend_inline.set_matplotlib_formats('svg')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Initialization" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We would like to couple MD particles that interact as nearly-hard spheres to MPCD particles in order to introduce approximate hydrodynamic interactions between them. HOOMD simulations are often initialized from particles saved in GSD files. MPCD particles cannot currently be added to a GSD file, but you can still initialize an MPCD simulation from a GSD file in two steps. We will illustrate how to do that here!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Making initial configuration for MD particles" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will start by making a configuration of MD particles that we write to a GSD file (Tutorial: [Initializing the System State](../00-Introducing-HOOMD-blue/03-Initializing-the-System-State.ipynb)). First, we place $N = \\rho V$ solute particles randomly on a lattice in a `Frame` box." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "solute_mass = 5\n", + "solute_density = 0.2\n", + "mpcd_density = 5\n", + "L = 20\n", + "kT = 1.0\n", + "\n", + "frame = gsd.hoomd.Frame()\n", + "rng = numpy.random.default_rng(seed=42)\n", + "\n", + "frame.configuration.box = [L, L, L, 0, 0, 0]\n", + "\n", + "frame.particles.N = numpy.round(solute_density * L**3).astype(int)\n", + "frame.particles.types = [\"A\"]\n", + "frame.particles.typeid = numpy.zeros(frame.particles.N, dtype=int)\n", + "solute_spacing = 1.2\n", + "cells = numpy.linspace(\n", + " -L / 2, L / 2, numpy.floor(L / solute_spacing).astype(int), endpoint=False\n", + ")\n", + "lattice = numpy.array(list(itertools.product(cells, repeat=3)))\n", + "chosen_lattice_sites = rng.choice(\n", + " a=lattice.shape[0], size=frame.particles.N, replace=False\n", + ")\n", + "frame.particles.position = numpy.take(lattice, chosen_lattice_sites, axis=0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "MD particles are usually heavier than the MPCD particles, and here we're setting them to have mass equal to the mass of MPCD particles in a collision cell because they are roughly the same size." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "frame.particles.mass = numpy.full(frame.particles.N, solute_mass)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We randomly draw velocities to be consistent with the Maxwell-Boltzmann distribution\n", + "for $k_{\\rm B} T = 1$ and their mass, and we force them to have zero mean. We save this `Frame` to a GSD file `init.gsd`. " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "vel = rng.normal(\n", + " loc=0.0, scale=numpy.sqrt(1.0 / mpcd_density), size=(frame.particles.N, 3)\n", + ")\n", + "vel -= numpy.mean(vel, axis=0)\n", + "frame.particles.velocity = vel\n", + "\n", + "with gsd.hoomd.open(name='init.gsd', mode='w') as f:\n", + " f.append(frame)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Adding MPCD particles" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will now add the MPCD particles using the `Snapshot` interface. We create `Snapshot` from the GSD file `init.gsd`." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "with gsd.hoomd.open('init.gsd') as traj:\n", + " snapshot = hoomd.Snapshot.from_gsd_frame(\n", + " traj[0], hoomd.communicator.Communicator()\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then, we will add the MPCD particles to the box at a number density of\n", + "$\\rho = 5$. This will be a bulk simulation, so the MPCD particles\n", + "should fill the box volume *V*. Hence, the number of MPCD particles to add\n", + "is $N = \\rho V$. We do this at random as in the\n", + "[pressure-driven flow](01-Pressure-Driven-Flow.ipynb) tutorial, assuming an\n", + "orthorhombic box." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "box = hoomd.Box.from_box(snapshot.configuration.box)\n", + "rng = numpy.random.default_rng(seed=42)\n", + "\n", + "snapshot.mpcd.types = [\"A\"]\n", + "snapshot.mpcd.N = numpy.round(mpcd_density * box.volume).astype(int)\n", + "snapshot.mpcd.position[:] = box.L * rng.uniform(low=-0.5,\n", + " high=0.5,\n", + " size=(snapshot.mpcd.N, 3)\n", + " )\n", + "\n", + "vel = rng.normal(loc=0.0, scale=numpy.sqrt(kT), size=(snapshot.mpcd.N, 3))\n", + "vel -= numpy.mean(vel, axis=0)\n", + "snapshot.mpcd.velocity[:] = vel" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The `Snapshot` is now fully initialized, so we can use it to create a\n", + "`Simulation`." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "sim = hoomd.Simulation(device=hoomd.device.CPU(), seed=1)\n", + "sim.create_state_from_snapshot(snapshot)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Configuring the MPCD integrator" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The MPCD `Integrator` handles the integration of both the MPCD particles (via\n", + "`collision_method` and `streaming_method`) and the MD particles (via `methods`\n", + "and `forces`). The MD particles are always integrated every timestep, but the\n", + "MPCD particles are only integrated on timesteps that are multiples of the\n", + "periods for streaming and collision.\n", + "\n", + "Hence, we first initialize the MPCD `Integrator` using the timestep that we\n", + "would normally use for an MD simulation because this timestep is shorter than\n", + "the MPCD streaming and collision times." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "integrator = hoomd.mpcd.Integrator(dt=0.005)\n", + "sim.operations.integrator = integrator" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will now perform the usual setup for an MD simulation, then configure\n", + "an MPCD simulation using the SRD collision method." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### MD integration method" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will define the interactions between MD particles as the repulsive Weeks-Chandler-Andersen (WCA) potential. These are attached to the integrator in the normal way." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "cell = hoomd.md.nlist.Cell(buffer=0.4)\n", + "wca = hoomd.md.pair.LJ(nlist=cell)\n", + "wca.params[(\"A\", \"A\")] = dict(epsilon=1, sigma=1)\n", + "wca.r_cut[(\"A\", \"A\")] = 2 ** (1.0 / 6.0)\n", + "integrator.forces.append(wca)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will use a `ConstantVolume` (NVE) integration method for the MD particles.\n", + "This integration method is usually recommended for MD particles\n", + "that are coupled to MPCD particles because the MPCD particles will act like a\n", + "temperature bath." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "nve = hoomd.md.methods.ConstantVolume(filter=hoomd.filter.All())\n", + "integrator.methods.append(nve)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### MPCD methods" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will now setup the SRD collision method for the MPCD particles using\n", + "collision time $\\Delta t = 0.1$, collision angle $\\alpha = 130^\\circ$, and a\n", + "thermostat to maintain constant temperature, which is the same setup\n", + "as we used in the [pressure-driven flow tutorial](01-Pressure-Driven-Flow.ipynb).\n", + "\n", + "The main difference for the collision method is that the `period` is now longer\n", + "because there are 20 MD timesteps per collision time. The solvent sorting period\n", + "should be similar extended, and here we use every 20 collisions." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "integrator.collision_method = hoomd.mpcd.collide.StochasticRotationDynamics(\n", + " period=20, angle=130, kT=1.0, embedded_particles=hoomd.filter.All()\n", + ")\n", + "integrator.mpcd_particle_sorter = hoomd.mpcd.tune.ParticleSorter(\n", + " trigger=integrator.collision_method.period * 20\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Apply the `Bulk` streaming method to propogate solvent particles with no confining geometry.\n", + "Streaming is only performed every collision step because the particles move\n", + "with constant velocity between collisions. Note that this means that\n", + "the MD particles and MPCD particles will appear to be \"out of sync\" if a\n", + "`Snapshot` is taken at a timestep that is not a multiple of the streaming period." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "integrator.streaming_method = hoomd.mpcd.stream.Bulk(\n", + " period=integrator.collision_method.period\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run simulation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that the simulation is configured, we will we run a short simulation to\n", + "write the MD particles trajectories to a GSD file. We will save the MD particles\n", + "every 2000 steps." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "gsd_writer = hoomd.write.GSD(\n", + " trigger=2000,\n", + " filename=\"monomers.gsd\",\n", + " dynamic=[\"particles/position\", \"particles/image\"],\n", + ")\n", + "sim.operations.writers.append(gsd_writer)\n", + "sim.run(200000)\n", + "gsd_writer.flush()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Analysis\n", + "\n", + "This trajectory can be used to compute the mean squared displacement of the MD\n", + "particles. First we read in the trajectory and unwrap the particle positions using freud." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "with gsd.hoomd.open(\"monomers.gsd\") as traj:\n", + " num_frames = len(traj)\n", + " N = traj[0].particles.N\n", + " positions = numpy.zeros((num_frames, N, 3), dtype=float)\n", + " for i, snap in enumerate(traj):\n", + " box = freud.box.Box.from_box(snap.configuration.box)\n", + " positions[i] = box.unwrap(snap.particles.position, snap.particles.image)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we calculate and plot the mean squared displacement using freud." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " 2024-06-18T11:12:01.659809\n", + " image/svg+xml\n", + " \n", + " \n", + " Matplotlib v3.9.0, https://matplotlib.org/\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "\n" + ], + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "msd = freud.msd.MSD(box, mode=\"window\")\n", + "msd.compute(positions)\n", + "msd.plot();" + ] + } + ], + "metadata": { + "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.1.undefined" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/AUTHORS.md b/AUTHORS.md index 500b4dd..5f5e126 100644 --- a/AUTHORS.md +++ b/AUTHORS.md @@ -6,3 +6,6 @@ * Tobias Dwyer - University of Michigan * Brandon Butler - University of Michigan * Tim Moore - University of Michigan +* Jinny Cha - Auburn University +* Michael P. Howard - Auburn University +* C. Levi Petix - Auburn University diff --git a/README.md b/README.md index fa80527..b9497cd 100644 --- a/README.md +++ b/README.md @@ -17,6 +17,7 @@ included in HOOMD-blue's [documentation]. * [Organizing and Executing Simulations](05-Organizing-and-Executing-Simulations/00-index.ipynb) * [Modelling Rigid Bodies](06-Modelling-Rigid-Bodies/00-index.ipynb) * [Modelling Patchy Particles](07-Modelling-Patchy-Particles/00-index.ipynb) +* [Multiparticle Collision Dynamics](09-Multiparticle-Collision-Dynamics/00-index.ipynb) ## Executing the tutorials