diff --git a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb index e591dfb4ef2..9501412c6ff 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb @@ -9,11 +9,12 @@ "#### Before you start:\n", "\n", "With this tutorial you can get started using the lattice-Boltzmann method\n", - "for scientific applications. We give a brief introduction about the theory\n", + "for hydrodynamics. We give a brief introduction about the theory\n", "and how to use the method in **ESPResSo**. We have selected two interesting problems for\n", "which LB can be applied and which are well understood. You can start with any of them.\n", "\n", - "The tutorial is relatively long and working through it carefully is work for at least a full day. You can however get a glimpse of different aspects by starting to work on the tasks.\n", + "The tutorial is relatively long and working through it carefully will take at least a full day.\n", + "You can however get a glimpse of different aspects by starting to work on the tasks.\n", "\n", "Note: LB cannot be used as a black box. It is unavoidable to spend time\n", "learning the theory and gaining practical experience." @@ -28,14 +29,9 @@ "In this tutorial, you will learn the basics of the lattice-Boltzmann method (LBM) with\n", "special focus on the application on soft matter simulations or more precisely on how to\n", "apply it in combination with molecular dynamics to take into account hydrodynamic\n", - "solvent effects without the need to introduce thousands of solvent particles.\n", - "\n", - "The LBM – its theory as well as its applications – is still a very active field of research.\n", - "After almost 20 years of development there are many cases in which the LBM has proven\n", - "to be fruitful, in other cases the LBM is considered promising, and in some cases it has\n", - "not been of any help. We encourage you to contribute to the scientific discussion of\n", - "the LBM because there is still a lot that is unknown or only vaguely known about this\n", - "fascinating method.\n", + "solvent effects without the need to introduce thousands of solvent particles. The LBM\n", + "– its theory as well as its applications – is a very active field of research with\n", + "more than 30 years of development.\n", "\n", "### Tutorial Outline\n", "\n", @@ -46,7 +42,7 @@ "\n", "* **Polymer diffusion.** We analyze the length dependence of the diffusion of polymers.\n", "\n", - "* **Poiseuille flow.** We reproduce the flow profile between two walls.\n", + "* **Poiseuille flow.** We reproduce the flow profile of a fluid moving through a pipe under a homogeneous force density.\n", "\n", "### Notes on the **ESPResSo** version you will need\n", "\n", @@ -55,6 +51,7 @@ "\n", "For the tutorial you will have to compile in the following features:\n", "```c++\n", + "#define LB_BOUNDARIES\n", "#define LB_BOUNDARIES_GPU\n", "#define LENNARD_JONES\n", "```\n", @@ -96,7 +93,7 @@ "\n", "The LBM discretizes the Boltzmann equation not only in real space (the lattice!) and\n", "time, but also the velocity space is discretized. A surprisingly small number of velocities,\n", - "in 3D usually 19, is sufficient to describe incompressible, viscous flow correctly. Mostly\n", + "usually 19 in 3D, is sufficient to describe incompressible, viscous flow correctly. Mostly\n", "we will refer to the three-dimensional model with a discrete set of 19 velocities, which is\n", "conventionally called D3Q19. These velocities, $\\vec{c_i}$ , are chosen such that they correspond to\n", "the movement from one lattice node to another in one time step. A two step scheme is\n", @@ -149,7 +146,7 @@ "why the mode space is helpful. A 19 dimensional matrix that\n", "conserves the first 4 modes (with the eigenvalue 1) is diagonal in the\n", "first four rows and columns.\n", - "Some struggling with lattice symmetries shows that four independent\n", + "By symmetry consideration one finds that only four independent\n", "variables are enough to characterize the linear relaxation\n", "process so that all symmetries of the lattice are obeyed. \n", "Two of them are closely related to \n", @@ -171,9 +168,7 @@ "In other words: each mode is relaxed towards\n", "its equilibrium value with a relaxation rate $\\gamma_i$.\n", "The conserved modes are not relaxed, i.e. their relaxation rate is 1.\n", - "\n", - "By symmetry consideration one finds that only four independent\n", - "relaxation rates are allowed. We summarize them here.\n", + "We summarize them here:\n", "\n", "\\begin{align*}\n", " m^\\star_i &= \\gamma_i m_i \\\\\n", @@ -188,7 +183,8 @@ "random fluctuations are added to the nonconserved modes $5\\dots 19$ on every LB node so that\n", "the LB fluid temperature is well defined and the corresponding\n", "fluctuation formula holds, according to the fluctuation dissipation theorem.\n", - "An extensive discussion of this topic is found in [3]." + "An extensive discussion of this topic is found in [3].\n", + "Thermalization of the fluid is optional in **ESPResSo**." ] }, { @@ -197,7 +193,7 @@ "source": [ "### Particle coupling\n", "\n", - "In **ESPResSo** particles are coupled to the LB fluid via the so called the force coupling method:\n", + "In **ESPResSo** particles are coupled to the LB fluid via the so-called force coupling method:\n", "the fluid velocity $\\vec{u}$ at the position of a particle is calculated \n", "via trilinear interpolation and a force is applied on the particle\n", "that is proportional to the velocity difference between particle and fluid:\n", @@ -224,8 +220,12 @@ "source": [ "## 3 The LB interface in ESPResSo\n", "\n", - "**ESPResSo** features two virtually independent implementations of LB. One implementation uses CPUs and one uses a GPU to perform the computational work. For this, we provide two actor classes LBFluid and LBFluidGPU in the module espressomd.lb,\n", - "as well as the optional LBBoundary class found in espressomd.lbboundaries.\n", + "**ESPResSo** features two virtually independent implementations of LB. One implementation uses CPUs and one uses a GPU to perform the computational work. For this, we provide two actor classes\n", + "[LBFluid](http://espressomd.org/html/doc/espressomd.html#espressomd.lb.LBFluid) and\n", + "[LBFluidGPU](http://espressomd.org/html/doc/espressomd.html#espressomd.lb.LBFluidGPU) in the module\n", + "[espressomd.lb](http://espressomd.org/html/doc/espressomd.html#module-espressomd.lb), as well as the optional\n", + "[LBBoundary](http://espressomd.org/html/doc/espressomd.html#espressomd.lbboundaries.LBBoundary) class found in\n", + "[espressomd.lbboundaries](http://espressomd.org/html/doc/espressomd.html#module-espressomd.lbboundaries).\n", "\n", "The LB lattice is a cubic lattice, with a lattice constant agrid that\n", "is the same in all spatial directions. The chosen box length must be an integer multiple\n", @@ -268,9 +268,11 @@ "+ agrid: The lattice constant of the fluid. It is used to determine the number of LB nodes per direction from box_l. *They have to be compatible.*\n", "+ visc: The kinematic viscosity\n", "+ tau: The time step of LB. It has to be an integer multiple of System.time_step.\n", + "+ kT: Thermal energy of the simulated heat bath for thermalized fluids, use 0 to deactivate thermalization.\n", + "+ seed: The random number generator seed, only relevant for thermalized fluids (i.e. kT \\> 0).\n", "+ ext_force_density: An external force density applied to every node. This is given as a list, tuple or array with three components.\n", "\n", - "Using these arguments, one can initialize an LBFluid object. This object then needs to be added to the system's actor list. The code below provides a minimum example.\n", + "Using these arguments, one can initialize an LBFluid object. This object then needs to be added to the system's actor list. The code below provides a minimal example.\n", "\n", "```python\n", "import espressomd\n", @@ -281,7 +283,7 @@ "system.time_step = 0.01\n", "system.cell_system.skin = 0.4\n", "\n", - "# Initialize and LBFluid with the minimum set of valid parameters.\n", + "# Initialize an LBFluid with the minimum set of valid parameters.\n", "lbf = lb.LBFluidGPU(agrid=1, dens=10, visc=.1, tau=0.01)\n", "# Activate the LB by adding it to the System's actor list.\n", "system.actors.add(lbf)\n", @@ -313,15 +315,13 @@ "source": [ "### The LBBoundary class\n", "\n", - "The LBBoundary class represents a boundary on the LBFluid\n", - "lattice. It depends on the classes of the module espressomd.shapes as it derives its geometry from them. For the initialization, the arguments shape and velocity are supported. The shape argument takes an object from the shapes module and the velocity argument expects a list, tuple or array containing 3 floats. Setting the velocity will results in a slip boundary condition.\n", + "The [LBBoundary](http://espressomd.org/html/doc/espressomd.html#espressomd.lbboundaries.LBBoundary) class represents a boundary on the\n", + "[LBFluid](http://espressomd.org/html/doc/espressomd.html#espressomd.lb.LBFluid) lattice.\n", + "It depends on the classes of the module espressomd.shapes as it derives its geometry from them. For the initialization, the arguments shape and velocity are supported. The shape argument takes an object from the shapes module and the velocity argument expects a list, tuple or array containing 3 floats. Setting the velocity will result in a slip boundary condition.\n", "\n", "Note that the boundaries are not constructed through the periodic boundary. If, for example, one would set a sphere with its center in one of the corner of the boxes, a sphere fragment will be generated. To avoid this, make sure the sphere, or any other boundary, fits inside the central box.\n", "\n", - "This part of the LB implementation is still experimental, so please tell us about your experience with it. In general even the simple case of no-slip\n", - "boundary is still an important research topic in the LB community and in combination with point particle coupling not much experience exists. This means: do research on that topic, play around with parameters and figure out what happens.\n", - "\n", - "Boundaries are initialized by passing a shape object to the LBBoundary class. One way to initialize a wall is:\n", + "Boundaries are instantiated by passing a shape object to the LBBoundary class. Here is one way to construct a wall and add it to an existing `system` instance:\n", "\n", "```python\n", "import espressomd.lbboundaries\n", @@ -332,7 +332,7 @@ "system.lbboundaries.add(wall)\n", "```\n", "\n", - "Note that all used variables are inherited from previous examples. This will create a wall with a surface normal of $(1, 0, 0)$ at a distance of 1 from the origin of the coordinate system in direction of the normal vector. The wall exhibits a slip boundary condition with a velocity of $(0, 0, 0.01)$. For a no-slip boundary condition, leave out the velocity argument or set it to zero. Please refer to the user guide for a complete list of constraints.\n", + "This will create a wall with a surface normal of $(1, 0, 0)$ at a distance of 1 from the origin of the coordinate system in direction of the normal vector. The wall exhibits a slip boundary condition with a velocity of $(0, 0, 0.01)$. For a no-slip boundary condition, leave out the velocity argument or set it to zero. Please refer to the user guide for a complete list of constraints.\n", "\n", "In **ESPResSo** the so-called *link bounce back* method is implemented, where the effective hydrodynamic boundary is located midway between boundary and fluid node." ] diff --git a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb index 345b19bc626..d23a2e6484f 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb @@ -11,15 +11,20 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## 5 Polymer Diffusion\n", + "## Diffusion of a single particle\n", "\n", - "In these exercises we want to use the LBM-MD-Hybrid to reproduce a classic result of polymer physics: the dependence of the diffusion coefficient of a polymer on its chain length. If no hydrodynamic interactions are present, one expects a scaling law $D \\propto N ^{- 1}$ and if they are present, a scaling law $D \\propto N^{- \\nu}$ is expected. Here $\\nu$ is the Flory exponent\n", - "that plays a very prominent role in polymer physics. It has a value of $\\sim 3/5$ in good\n", - "solvent conditions in 3D. Discussions on these scaling laws can be found in polymer\n", - "physics textbooks like [4–6].\n", + "In these exercises we want to reproduce a classic result of polymer physics: the dependence \n", + "of the diffusion coefficient of a polymer on its chain length. If no hydrodynamic interactions\n", + "are present, one expects a scaling law $D \\propto N ^{- 1}$ and if they are present, a scaling law\n", + "$D \\propto N^{- \\nu}$ is expected. Here $\\nu$ is the Flory exponent that plays a very prominent\n", + "role in polymer physics. It has a value of $\\sim 3/5$ in good solvent conditions in 3D.\n", + "Discussions on these scaling laws can be found in polymer physics textbooks like [4–6].\n", "\n", - "The reason for the different scaling law is the following: when being transported, every monomer creates a flow field that follows the direction of its motion. This flow field makes it easier for other monomers to follow its motion. This makes a polymer (given it is sufficiently long) diffuse more like a compact object including the fluid inside it, although it does not have clear boundaries. It can be shown that its motion can be described by its\n", - "hydrodynamic radius. It is defined as:\n", + "The reason for the different scaling law is the following: when being transported, every monomer\n", + "creates a flow field that follows the direction of its motion. This flow field makes it easier for\n", + "other monomers to follow its motion. This makes a polymer (given it is sufficiently long) diffuse\n", + "more like a compact object including the fluid inside it, although it does not have clear boundaries.\n", + "It can be shown that its motion can be described by its hydrodynamic radius. It is defined as:\n", "\n", "\\begin{equation}\n", " \\left\\langle \\frac{1}{R_h} \\right\\rangle = \\left\\langle \\frac{1}{N^2}\\sum_{i\\neq j} \\frac{1}{\\left| r_i - r_j \\right|} \\right\\rangle\n", @@ -29,7 +34,7 @@ "and the diffusion coefficient of a long polymer is proportional to its inverse $R_h$.\n", "\n", "For shorter polymers there is a transition region. It can be described\n", - "by the Kirkwood-Zimm model:\n", + "by the Kirkwood–Zimm model:\n", "\n", "\\begin{equation}\n", " D=\\frac{D_0}{N} + \\frac{k_B T}{6 \\pi \\eta } \\left\\langle \\frac{1}{R_h} \\right\\rangle\n", @@ -37,7 +42,7 @@ "\n", "Here $D_0$ is the monomer diffusion coefficient and $\\eta$ the \n", "viscosity of the fluid. For a finite system size the second part of the\n", - "diffusion is subject of a $1/L$ finite size effect, because\n", + "diffusion is subject to a $1/L$ finite size effect, because\n", "hydrodynamic interactions are proportional to the inverse\n", "distance and thus long ranged. It can be taken into account\n", "by a correction:\n", @@ -51,7 +56,7 @@ "too disappointed if you don't manage to do so.\n", "\n", "We want to determine the long-time self diffusion coefficient from the mean square\n", - "displacement of the center-of-mass a single polymer. For large $t$ the mean square displacement is\n", + "displacement of the center-of-mass of a single polymer. For large $t$ the mean square displacement is\n", "proportional to the time and the diffusion coefficient occurs as a \n", "prefactor:\n", "\n", @@ -59,9 +64,8 @@ " D = \\lim_{t\\to\\infty}\\left[ \\frac{1}{6t} \\left\\langle \\left(\\vec{r}(t) - \\vec{r}(0)\\right)^2 \\right\\rangle \\right].\n", "\\end{equation}\n", "\n", - "This equation can be\n", - "found in virtually any simulation textbook, like [7]. We will therefore set up a polymer\n", - "in an LB fluid, simulate for an appropriate amount of time, calculate the mean square\n", + "This equation can be found in virtually any simulation textbook, like [7]. We will set up a\n", + "polymer in an implicit solvent, simulate for an appropriate amount of time, calculate the mean square\n", "displacement as a function of time and obtain the diffusion coefficient from a linear\n", "fit. However we will have a couple of steps inbetween and divide the full problem into\n", "subproblems that allow to (hopefully) fully understand the process." @@ -71,14 +75,55 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### 5.1 Step 1: Diffusion of a single particle" + "### 1. Setting up the observable" ] }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "Write a function with signature `correlator_msd(pid, tau_max)` that returns a\n", + "mean-squared displacement correlator that is updated every time step." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "def correlator_msd(pid, tau_max):\n", + " pos = espressomd.observables.ParticlePositions(ids=(pid,))\n", + " pos_cor = espressomd.accumulators.Correlator(\n", + " obs1=pos, tau_lin=16, tau_max=tau_max, delta_N=1,\n", + " corr_operation=\"square_distance_componentwise\", compress1=\"discard1\")\n", + " return pos_cor\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Our first step is to investigate the diffusion of a single particle that is coupled to an LB fluid by the point coupling method:" + "### 2. Simulating the Brownian motion" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We will simulate the diffusion of a single particle that is coupled to an implicit solvent." ] }, { @@ -93,55 +138,44 @@ "\n", "import espressomd\n", "import espressomd.accumulators\n", - "import espressomd.lb\n", "import espressomd.observables\n", "\n", "logging.basicConfig(level=logging.INFO, stream=sys.stdout)\n", "\n", - "\n", "# Constants\n", - "LOOPS = 40000\n", - "STEPS = 10\n", + "KT = 1.1\n", + "STEPS = 400000\n", "\n", "# System setup\n", "system = espressomd.System(box_l=[16] * 3)\n", "system.time_step = 0.01\n", "system.cell_system.skin = 0.4\n", "\n", - "\n", - "lbf = espressomd.lb.LBFluidGPU(kT=1, seed=132, agrid=1, dens=1, visc=5, tau=0.01)\n", - "system.actors.add(lbf)\n", - "\n", - "\n", "system.part.add(pos=[0, 0, 0])\n", "\n", - "# Setup observable\n", - "pos = espressomd.observables.ParticlePositions(ids=(0,))\n", - "\n", - "# Run for different friction constants\n", - "lb_gammas = [1.0, 2.0, 4.0, 10.0]\n", + "# Run for different friction coefficients\n", + "gammas = [1.0, 2.0, 4.0, 10.0]\n", + "tau_results = []\n", "msd_results = []\n", "\n", - "for gamma in lb_gammas:\n", + "for gamma in gammas:\n", " system.auto_update_accumulators.clear()\n", " system.thermostat.turn_off()\n", - " system.thermostat.set_lb(LB_fluid=lbf, seed=123, gamma=gamma)\n", - " # Equilibrate\n", + " system.thermostat.set_langevin(kT=KT, gamma=gamma, seed=42)\n", + "\n", " logging.info(\"Equilibrating the system.\")\n", " system.integrator.run(1000)\n", " logging.info(\"Equilibration finished.\")\n", + "\n", " # Setup observable correlator\n", - " correlator = espressomd.accumulators.Correlator(obs1=pos, tau_lin=16, tau_max=LOOPS * STEPS, delta_N=1,\n", - " corr_operation=\"square_distance_componentwise\", compress1=\"discard1\")\n", + " correlator = correlator_msd(0, STEPS)\n", " system.auto_update_accumulators.add(correlator)\n", "\n", " logging.info(\"Sampling started for gamma = {}.\".format(gamma))\n", - " for i in range(LOOPS):\n", - " system.integrator.run(STEPS)\n", + " system.integrator.run(STEPS)\n", " correlator.finalize()\n", - " msd_results.append(np.column_stack((correlator.lag_times(),\n", - " correlator.sample_sizes(),\n", - " correlator.result().reshape([-1, 3]))))\n", + " tau_results.append(correlator.lag_times())\n", + " msd_results.append(np.sum(correlator.result().reshape([-1, 3]), axis=1))\n", "\n", "logging.info(\"Sampling finished.\")" ] @@ -150,7 +184,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### Plotting the results" + "### 3. Data analysis\n", + "#### 3.1 Plotting the results" ] }, { @@ -164,16 +199,14 @@ "\n", "plt.rcParams.update({'font.size': 22})\n", "\n", - "plt.figure(figsize=(10,10))\n", - "plt.xlabel('tau [$\\Delta t$]')\n", - "plt.ylabel('msd [$\\sigma^2$]')\n", - "for index, r in enumerate(msd_results):\n", - " # adding up componentwise MSDs\n", + "plt.figure(figsize=(10, 10))\n", + "plt.xlabel(r'$\\tau$ [$\\Delta t$]')\n", + "plt.ylabel(r'MSD [$\\sigma^2$]')\n", + "for index, (tau, msd) in enumerate(zip(tau_results, msd_results)):\n", " # We skip the first entry since it's zero by definition and cannot be displayed\n", " # in a loglog plot. Furthermore, we only look at the first 100 entries due to\n", - " # lack of good statistics for larger times.\n", - " msd = r[1:100, 2] + r[1:100, 3] + r[1:100, 4]\n", - " plt.loglog(r[1:100, 0], msd, label=r\"$\\gamma=${}\".format(str(lb_gammas[index])))\n", + " # the high variance for larger lag times.\n", + " plt.loglog(tau[1:100], msd[1:100], label=r'$\\gamma=${:.1f}'.format(gammas[index]))\n", "plt.legend()\n", "plt.show()" ] @@ -182,12 +215,145 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### Calculating the diffusion coefficient\n", - "In this script an LB fluid and a single particle are created and thermalized. The\n", - "random forces on the particle and within the LB fluid will cause the particle to move.\n", + "#### 3.2 Calculating the diffusion coefficient\n", + "\n", + "In this script an implicit solvent and a single particle are created and thermalized.\n", + "The random forces on the particle will cause the particle to move.\n", "The mean squared displacement is calculated during the simulation via a multiple-tau\n", "correlator. \n", - "Can you give an explanation for the quadratic time dependency for short times? Use the function curve_fit from the module scipy.optimize to produce a fit for the linear regime and determine the diffusion coefficients for the different $\\gamma$s. Calculate the diffusion coefficient for all cases and plot them as a function of γ. What relation do you observe?" + "Can you give an explanation for the quadratic time dependency for short times?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The MSD of a Brownian motion can be decomposed in three main regimes [8]:\n", + "* for short lag times $\\tau < \\tau_p$, the particle motion is not\n", + " significantly impeded by solvent collisions: it's in the ballistic mode\n", + " (collision-free regime) where $\\operatorname{MSD}(t) \\sim (k_BT / m) t^2$\n", + "* for long lag times $\\tau > \\tau_f$, the particle motion is determined by\n", + " numerous collisions with the solvent: it's in the diffusive mode where\n", + " $\\operatorname{MSD}(t) \\sim 6t$\n", + "* for lag times between $\\tau_p$ and $\\tau_f$, there is a crossover mode\n", + "\n", + "The values $\\tau_p$ and $\\tau_f$ can be obtained manually through visual\n", + "inspection of the MSD plot, or more accurately by non-linear fitting [9].\n", + "\n", + "The cutoff lag time $\\tau_p$ between the ballistic and crossover modes is proportional\n", + "to the particle mass and inversely proportional to the friction coefficient. \n", + "In the graph below, a parabola is fitted to the data points in the ballistic mode for\n", + "each $\\gamma$ and plotted beyond the crossover region to reveal the deviation from the\n", + "ballistic mode. This deviation is clearly visible in the $\\gamma = 10$ case, because\n", + "the assumption of a collision-free regime quickly breaks down when a particle is\n", + "coupled to its surrounding fluid with a high friction coefficient." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import scipy.optimize\n", + "\n", + "def quadratic(x, a, b, c):\n", + " return a * x**2 + b * x + c\n", + "\n", + "# cutoffs for the ballistic regime (different for each gamma value)\n", + "tau_p_values = [14, 12, 10, 7]\n", + "\n", + "plt.figure(figsize=(10, 10))\n", + "plt.xlabel(r'$\\tau$ [$\\Delta t$]')\n", + "plt.ylabel(r'MSD [$\\sigma^2$]')\n", + "for index, (tau_p, tau, msd) in enumerate(zip(tau_p_values, tau_results, msd_results)):\n", + " (a, b, c), _ = scipy.optimize.curve_fit(quadratic, tau[:tau_p], msd[:tau_p])\n", + " x = np.linspace(tau[0], tau[max(tau_p_values) - 1], 50)\n", + " p = plt.plot(x, quadratic(x, a, b, c), '-')\n", + " plt.plot(tau[:max(tau_p_values)], msd[:max(tau_p_values)], 'o', color=p[0].get_color(),\n", + " label=r'$\\gamma=${:.1f}'.format(gammas[index]))\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Use the function [curve_fit()](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) from the module scipy.optimize to produce a fit for the linear regime and determine the diffusion coefficients for the different $\\gamma$s." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For large $t$ the diffusion coefficient can be expressed as:\n", + "\n", + "$$6D = \\lim_{t\\to\\infty} \\frac{\\partial \\operatorname{MSD}(t)}{\\partial t}$$\n", + "\n", + "which is simply the slope of the MSD in the diffusive mode." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def linear(x, a, b):\n", + " return a * x + b\n", + "\n", + "# cutoffs for the diffusive regime (different for each gamma value)\n", + "tau_f_values = [24, 22, 20, 17]\n", + "# cutoff for the data series (larger lag times have larger variance due to undersampling)\n", + "cutoff_limit = 90\n", + "\n", + "diffusion_results = []\n", + "\n", + "plt.figure(figsize=(10, 8))\n", + "plt.xlabel(r'$\\tau$ [$\\Delta t$]')\n", + "plt.ylabel(r'MSD [$\\sigma^2$]')\n", + "for index, (tau_f, tau, msd) in enumerate(zip(tau_f_values, tau_results, msd_results)):\n", + " (a, b), _ = scipy.optimize.curve_fit(linear, tau[tau_f:cutoff_limit], msd[tau_f:cutoff_limit])\n", + " x = np.linspace(tau[tau_f], tau[cutoff_limit - 1], 50)\n", + " p = plt.plot(x, linear(x, a, b), '-')\n", + " plt.plot(tau[tau_f:cutoff_limit], msd[tau_f:cutoff_limit], 'o', color=p[0].get_color(),\n", + " label=r'$\\gamma=${:.1f}'.format(gammas[index]))\n", + " diffusion_results.append(a / 6)\n", + "\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Calculate the diffusion coefficient for all cases and plot them as a function of $\\gamma$. What relation do you observe?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the diffusive mode, one can derive $D = k_BT / \\gamma$ from the Stokes–Einstein relation [8]." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(10, 8))\n", + "plt.xlabel(r'$\\gamma$')\n", + "plt.ylabel('Diffusion coefficient [$\\sigma^2/t$]')\n", + "x = np.linspace(0.9 * min(gammas), 1.1 * max(gammas), 50)\n", + "y = KT / x\n", + "plt.plot(x, y, '-', label=r'$k_BT\\gamma^{-1}$')\n", + "plt.plot(gammas, diffusion_results, 'o', label='D')\n", + "plt.legend()\n", + "plt.show()" ] }, { @@ -202,7 +368,9 @@ "[4] P. G. de Gennes. *Scaling Concepts in Polymer Physics*. Cornell University Press, Ithaca, NY, 1979. \n", "[5] M. Doi. *Introduction to Polymer Physics.* Clarendon Press, Oxford, 1996. \n", "[6] Michael Rubinstein and Ralph H. Colby. *Polymer Physics.* Oxford University Press, Oxford, UK, 2003. \n", - "[7] Daan Frenkel and Berend Smit. *Understanding Molecular Simulation.* Academic Press, San Diego, second edition, 2002." + "[7] Daan Frenkel and Berend Smit. *Understanding Molecular Simulation*, section 4.4.1. Academic Press, San Diego, second edition, 2002. \n", + "[8] R. Huang, I. Chavez, K. Taute, et al. Direct observation of the full transition from ballistic to diffusive Brownian motion in a liquid. *Nature Phys.*, 7, 2011. doi:[10.1038/nphys1953](https://doi.org/10.1038/nphys1953) \n", + "[9] M. K. Riahi, I. A. Qattan, J. Hassan, D. Homouz, Identifying short- and long-time modes of the mean-square displacement: An improved nonlinear fitting approach. *AIP Advances*, 9:055112, 2019. doi:[10.1063/1.5098051](https://doi.org/10.1063/1.5098051) " ] } ], diff --git a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb index 93a7986821c..12a76ef4074 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb @@ -11,23 +11,21 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### 5.2 Step 2: Diffusion of a polymer\n", + "## Diffusion of a polymer\n", "\n", "One of the typical applications of **ESPResSo** is the simulation of polymer chains with a bead-spring-model. For this we need a repulsive interaction between all beads, for which one usually takes a shifted and truncated Lennard-Jones (so-called Weeks–Chandler–Andersen or WCA) interaction, and additionally a bonded interaction between adjacent beads to hold the polymer together. You have already learned that the command\n", "\n", "```python\n", "system.non_bonded_inter[0, 0].lennard_jones.set_params(\n", - " epsilon=1.0, sigma=1.0,\n", - " shift=0.25, cutof=1.226)\n", + " epsilon=1.0, sigma=1.0, shift=0.25, cutoff=1.1225)\n", "```\n", "\n", "creates a Lennard-Jones interaction with $\\varepsilon=1.$, $\\sigma=1.$,\n", - "$r_\\text{cut} = 1.125$ and $\\varepsilon_\\text{shift}=0.25$ between particles\n", - "of type 0, which is the desired \n", - "repulsive interaction. The command\n", + "$r_\\text{cut} = 1.1225$ and $\\varepsilon_\\text{shift}=0.25$ between particles\n", + "of type 0, which is the desired repulsive interaction. The command\n", "\n", "```python\n", - "fene = espressomd.interactions.FeneBond(k=7, d_r_max=2)\n", + "fene = espressomd.interactions.FeneBond(k=7, r_0=1, d_r_max=2)\n", "```\n", "\n", "creates a FeneBond object (see **ESPResSo** manual for the details). What is left to be done is to add this bonded interaction to the system via\n", @@ -44,15 +42,148 @@ "```python\n", "positions = espressomd.polymer.linear_polymer_positions(n_polymers=1,\n", " beads_per_chain=10,\n", - " bond_length=1, seed=5642,\n", + " bond_length=1, seed=42,\n", " min_distance=0.9)\n", "```\n", "\n", - "which would create positions for a single polymer with 10 monomers. Please check the documentation for a more detailed description.\n", + "which would create positions for a single polymer with 10 monomers. Please check the documentation for a more detailed description." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1. Setting up the polymer and observables\n", + "\n", + "The first task is to compute the average hydrodynamic radius $R_h$, end-to-end distance $R_F$\n", + "and radius of gyration $R_g$ for different polymer lengths. This will be achieved with the\n", + "corresponding observables described in the user guide under\n", + "[Analysis / Direct analysis routines / Chains](http://espressomd.org/html/doc/analysis.html#chains).\n", "\n", - "Furthermore we want to compute the diffusion constant of the polymer for different numbers of monomers. For this purpose we can again use the multiple tau correlator. Have a look at the **ESPResSo**-script for the single particle diffusion. Find out how many integration steps are necessary to capture the long-time diffusion regime of the polymer. The following script computes the mean squared displacement for the center of mass of the polymer as well as the average hydrodynamic radius. \n", + "The second task is to estimate the polymer diffusion coefficient for different polymer lengths\n", + "using two methods:\n", + "* the center of mass mean squared displacement method (introduced in a previous part of this tutorial)\n", + "* the center of mass velocity autocorrelation method (also known as Green–Kubo method)\n", "\n", - "Adapt the script for sampling different numbers of monomers and analyze the result. How does the hydrodynamic radius and the diffusion constant depend on the number of monomers? A solution can be found in doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part3_solution.py." + "For this purpose we can again use the [multiple tau correlator](http://espressomd.org/html/doc/analysis.html#details-of-the-multiple-tau-correlation-algorithm)." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "Write a function with signature `build_polymer(system, n_monomers, polymer_params, fene)` that creates\n", + "a linear polymer made of `n_monomers` particles, with parameters `polymer_params`. The particles need\n", + "to be created and linked together with the `fene` bond." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "def build_polymer(system, n_monomers, polymer_params, fene):\n", + " positions = espressomd.polymer.linear_polymer_positions(\n", + " beads_per_chain=n_monomers, **polymer_params)\n", + " for i, pos in enumerate(positions[0]):\n", + " pid = len(system.part)\n", + " system.part.add(id=pid, pos=pos)\n", + " if i > 0:\n", + " system.part[pid].add_bond((fene, pid - 1))\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "Write a function with signature `correlator_msd(pids_monomers, tau_max)` that returns a center-of-mass\n", + "mean-squared displacement correlator that is updated every time step, and a function with signature\n", + "`correlator_gk(pids_monomers, tau_max)` that returns a center-of-mass velocity correlator that is updated\n", + "every 10 time steps. You can find exemples in the user guide section\n", + "[calculating a particle's diffusion coefficient](http://espressomd.org/html/doc/analysis.html#example-calculating-a-particle-s-diffusion-coefficient)." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "def correlator_msd(pids_monomers, tau_max):\n", + " com_pos = espressomd.observables.ComPosition(ids=pids_monomers)\n", + " com_pos_cor = espressomd.accumulators.Correlator(\n", + " obs1=com_pos, tau_lin=16, tau_max=tau_max, delta_N=1,\n", + " corr_operation=\"square_distance_componentwise\", compress1=\"discard1\")\n", + " return com_pos_cor\n", + "\n", + "def correlator_gk(pids_monomers, tau_max):\n", + " com_vel = espressomd.observables.ComVelocity(ids=pids_monomers)\n", + " com_vel_cor = espressomd.accumulators.Correlator(\n", + " obs1=com_vel, tau_lin=16, tau_max=tau_max, delta_N=10,\n", + " corr_operation=\"scalar_product\", compress1=\"discard1\")\n", + " return com_vel_cor\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can simulate a polymer in the Rouse regime using an implicit solvent model, e.g. Langevin dynamics,\n", + "or in the Zimm regime using a lattice-Boltzmann fluid." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def solvent_langevin(system, kT, gamma):\n", + " '''\n", + " Implicit solvation model based on Langevin dynamics (Rouse model).\n", + " '''\n", + " system.thermostat.set_langevin(kT=kT, gamma=gamma, seed=42)\n", + "\n", + "def solvent_lbm(system, kT, gamma):\n", + " '''\n", + " Lattice-based solvation model based on the LBM (Zimm model).\n", + " '''\n", + " lbf = espressomd.lb.LBFluidGPU(kT=kT, seed=42, agrid=1, dens=1,\n", + " visc=5, tau=system.time_step)\n", + " system.actors.add(lbf)\n", + " system.thermostat.set_lb(LB_fluid=lbf, gamma=gamma, seed=42)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2. Simulating the polymer" ] }, { @@ -75,83 +206,577 @@ "\n", "espressomd.assert_features(['LENNARD_JONES'])\n", "\n", - "# Setup constant\n", + "# Setup constants\n", "TIME_STEP = 0.01\n", - "LOOPS = 100\n", + "LOOPS = 4000\n", "STEPS = 100\n", + "KT = 1.0\n", + "GAMMA = 5.0\n", + "POLYMER_PARAMS = {'n_polymers': 1, 'bond_length': 1, 'seed': 42, 'min_distance': 0.9}\n", + "POLYMER_MODEL = 'Rouse'\n", + "assert POLYMER_MODEL in ('Rouse', 'Zimm')\n", + "if POLYMER_MODEL == 'Zimm':\n", + " espressomd.assert_features(['CUDA'])\n", + " import espressomd.lb\n", "\n", "# System setup\n", - "system = espressomd.System(box_l=[32.0, 32.0, 32.0])\n", + "system = espressomd.System(box_l=[12.0, 12.0, 12.0])\n", "system.cell_system.skin = 0.4\n", "\n", - "N_MONOMERS = 20 # The number of monomers has been set to be 20 as default\n", - " # Change this value for further simulations\n", - "\n", "# Lennard-Jones interaction\n", "system.non_bonded_inter[0, 0].lennard_jones.set_params(\n", - " epsilon=1.0, sigma=1.0,\n", - " shift=\"auto\", cutoff=2.0**(1.0 / 6.0))\n", + " epsilon=1.0, sigma=1.0, shift=\"auto\", cutoff=2.0**(1.0 / 6.0))\n", "\n", "# Fene interaction\n", - "fene = espressomd.interactions.FeneBond(k=7, d_r_max=2)\n", + "fene = espressomd.interactions.FeneBond(k=7, r_0=1, d_r_max=2)\n", "system.bonded_inter.add(fene)\n", "\n", + "N_MONOMERS = np.array([6, 8, 10])\n", + "\n", + "com_pos_tau_results = []\n", + "com_pos_msd_results = []\n", + "com_vel_tau_results = []\n", + "com_vel_acf_results = []\n", + "rh_results = []\n", + "rf_results = []\n", + "rg_results = []\n", + "for index, N in enumerate(N_MONOMERS):\n", + " logging.info(\"Polymer size: {}\".format(N))\n", + " build_polymer(system, N, POLYMER_PARAMS, fene)\n", + "\n", + " logging.info(\"Warming up the polymer chain.\")\n", + " system.time_step = 0.002\n", + " system.integrator.set_steepest_descent(\n", + " f_max=1.0,\n", + " gamma=10,\n", + " max_displacement=0.01)\n", + " system.integrator.run(2000)\n", + " system.integrator.set_vv()\n", + " logging.info(\"Warmup finished.\")\n", + "\n", + " logging.info(\"Equilibration.\")\n", + " system.time_step = TIME_STEP\n", + " system.thermostat.set_langevin(kT=1.0, gamma=50, seed=42)\n", + " system.integrator.run(2000)\n", + " logging.info(\"Equilibration finished.\")\n", + "\n", + " system.thermostat.turn_off()\n", + "\n", + " if POLYMER_MODEL == 'Rouse':\n", + " solvent_langevin(system, KT, GAMMA)\n", + " elif POLYMER_MODEL == 'Zimm':\n", + " solvent_lbm(system, KT, GAMMA)\n", + "\n", + " logging.info(\"Warming up the system with the fluid.\")\n", + " system.integrator.run(1000)\n", + " logging.info(\"Warming up the system with the fluid finished.\")\n", + "\n", + " # configure MSD correlator\n", + " com_pos_cor = correlator_msd(np.arange(N), LOOPS * STEPS)\n", + " system.auto_update_accumulators.add(com_pos_cor)\n", + "\n", + " # configure Green-Kubo correlator\n", + " com_vel_cor = correlator_gk(np.arange(N), LOOPS * STEPS)\n", + " system.auto_update_accumulators.add(com_vel_cor)\n", + "\n", + " logging.info(\"Sampling started.\")\n", + " rhs = np.zeros(LOOPS)\n", + " rfs = np.zeros(LOOPS)\n", + " rgs = np.zeros(LOOPS)\n", + " for i in range(LOOPS):\n", + " system.integrator.run(STEPS)\n", + " rhs[i] = system.analysis.calc_rh(\n", + " chain_start=0,\n", + " number_of_chains=1,\n", + " chain_length=N)[0]\n", + " rfs[i] = system.analysis.calc_re(\n", + " chain_start=0,\n", + " number_of_chains=1,\n", + " chain_length=N)[0]\n", + " rgs[i] = system.analysis.calc_rg(\n", + " chain_start=0,\n", + " number_of_chains=1,\n", + " chain_length=N)[0]\n", + " logging.info(\"Sampling finished.\")\n", + "\n", + " # store results\n", + " com_pos_cor.finalize()\n", + " com_pos_tau_results.append(com_pos_cor.lag_times())\n", + " com_pos_msd_results.append(np.sum(com_pos_cor.result(), axis=1))\n", + " com_vel_cor.finalize()\n", + " com_vel_tau_results.append(com_vel_cor.lag_times())\n", + " com_vel_acf_results.append(com_vel_cor.result())\n", + " rh_results.append(rhs)\n", + " rf_results.append(rfs)\n", + " rg_results.append(rgs)\n", + "\n", + " # reset system\n", + " system.part.clear()\n", + " system.thermostat.turn_off()\n", + " system.actors.clear()\n", + " system.auto_update_accumulators.clear()\n", + "\n", + "rh_results = np.array(rh_results)\n", + "rf_results = np.array(rf_results)\n", + "rg_results = np.array(rg_results)\n", + "com_pos_tau_results = np.array(com_pos_tau_results)\n", + "com_pos_msd_results = np.reshape(com_pos_msd_results, [len(N_MONOMERS),-1])\n", + "com_vel_tau_results = np.array(com_vel_tau_results)\n", + "com_vel_acf_results = np.reshape(com_vel_acf_results, [len(N_MONOMERS),-1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 3. Data analysis\n", "\n", - "# Setup polymer of part_id 0 with fene bond\n", - "positions = espressomd.polymer.linear_polymer_positions(n_polymers=1,\n", - " beads_per_chain=N_MONOMERS,\n", - " bond_length=1, seed=5642,\n", - " min_distance=0.9)\n", - "for i, pos in enumerate(positions[0]):\n", - " id = len(system.part)\n", - " system.part.add(id=id, pos=pos)\n", - " if i > 0:\n", - " system.part[id].add_bond((fene, id - 1))\n", + "We will calculate the means of time series with error bars obtained from\n", + "the correlation-corrected standard error of the mean [8,9]." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib notebook\n", + "import matplotlib.pyplot as plt\n", + "import matplotlib.ticker as ticker\n", + "plt.rcParams.update({'font.size': 22})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import scipy.signal\n", + "\n", + "def autocorrelation(time_series):\n", + " '''\n", + " Normalized autocorrelation function with zero-padding. This is a\n", + " scipy implementation of the algorithm described in de Buyl 2018,\n", + " JOSS 3(28), pp. 877 (https://doi.org/10.21105/joss.00877).\n", + " '''\n", + " n = time_series.size\n", + " n_with_padding = 2**int(np.ceil(np.log2(n)) + 1)\n", + " signal = np.zeros(n_with_padding)\n", + " signal[:n] = time_series\n", + " acf = scipy.signal.correlate(signal, signal, mode='full', method='fft')[-n_with_padding:][:n]\n", + " acf_normalized = acf / (n - np.arange(n)) # rescale by bin sizes\n", + " return acf_normalized\n", + "\n", + "def standard_error_mean_autocorrelation(time_series, variable_label):\n", + " '''\n", + " Calculate the mean and the correlation-corrected standard error\n", + " of the mean of time series by integrating the autocorrelation\n", + " function. See Janke 2002 [8] and Weigel, Janke 2010 [9].\n", + "\n", + " Due to the short simulation length, it is not possible to fit an\n", + " exponential to the long-time tail. Instead, return a percentile.\n", + " '''\n", + " summary = []\n", + " fig = plt.figure(figsize=(10, 6))\n", + " for signal, N in zip(time_series, N_MONOMERS):\n", + " acf = autocorrelation(signal - np.mean(signal))\n", + " # the acf cannot be integrated beyond tau=N/2\n", + " integral = np.array([acf[0] + 2 * np.sum(acf[1:j]) for j in np.arange(1, len(acf)//2)])\n", + " # remove the noisy part of the integral\n", + " negative_number_list = np.nonzero(integral < 0)\n", + " if negative_number_list[0].size:\n", + " integral = integral[:int(0.95 * negative_number_list[0][0])]\n", + " # compute the standard error of the mean\n", + " std_err = np.sqrt(integral / acf.size)\n", + " # due to the small sample size, the long-time tail is not\n", + " # well resolved and cannot be fitted, so we use a percentile\n", + " asymptote = np.percentile(std_err, 75)\n", + " # plot the integral and asymptote\n", + " p = plt.plot([0, len(std_err)], 2 * [asymptote], '--')\n", + " plt.plot(np.arange(len(std_err)) + 1, std_err,\n", + " '-', color=p[0].get_color(),\n", + " label=rf'$\\int {variable_label}$ for N={N}')\n", + " summary.append((np.mean(signal), asymptote))\n", + " plt.xlabel(r'Lag time $\\tau / \\Delta t$')\n", + " plt.ylabel(rf'$\\int_{{-\\tau}}^{{+\\tau}} {variable_label}$')\n", + " plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 0))\n", + " plt.legend()\n", + " plt.show()\n", + " return np.array(summary)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 3.1 Distance-based macromolecular properties\n", "\n", + "How do $R_h$, $R_g$, $R_F$ and the diffusion coefficient $D$ depend on the number of monomers?\n", + "You can refer to the Flory theory of polymers, and assume you are simulating a real polymer in a\n", + "good solvent, with Flory exponent $\\nu \\approx 0.588$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the end-to-end distance $R_F$ of the polymer as a function of the number of monomers. What relation do you observe?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The end-to-end distance follows the law $R_F = c_F N^\\nu$ with $c_F$ a constant and $\\nu$ the Flory exponent." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rf_summary = standard_error_mean_autocorrelation(rf_results, r'\\operatorname{acf}(R_F)')\n", + "rf_exponent, rf_prefactor = np.polyfit(np.log(N_MONOMERS), np.log(rf_summary[:,0]), 1)\n", + "rf_prefactor = np.exp(rf_prefactor)\n", + "\n", + "fig = plt.figure(figsize=(10, 6))\n", + "x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)\n", + "plt.plot(x, rf_prefactor * x**rf_exponent, '-',\n", + " label=rf'$R_F^{{\\mathrm{{fit}}}} = {rf_prefactor:.2f} N^{{{rf_exponent:.2f}}}$')\n", + "plt.errorbar(N_MONOMERS, rf_summary[:,0],\n", + " yerr=rf_summary[:,1],\n", + " ls='', marker='o', capsize=5, capthick=1,\n", + " label=r'$R_F^{\\mathrm{simulation}}$')\n", + "plt.xlabel('Number of monomers $N$')\n", + "plt.ylabel(r'End-to-end distance [$\\sigma$]')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the radius of gyration $R_g$ of the polymer as a function of the number of monomers. What relation do you observe?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The radius of gyration follows the law $R_g = c_g N^\\nu$ with $c_g$ a constant and $\\nu$ the Flory exponent." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rg_summary = standard_error_mean_autocorrelation(rg_results, r'\\operatorname{acf}(R_g)')\n", + "rg_exponent, rg_prefactor = np.polyfit(np.log(N_MONOMERS), np.log(rg_summary[:,0]), 1)\n", + "rg_prefactor = np.exp(rg_prefactor)\n", + "\n", + "fig = plt.figure(figsize=(10, 6))\n", + "x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)\n", + "plt.plot(x, rg_prefactor * x**rg_exponent, '-',\n", + " label=rf'$R_g^{{\\mathrm{{fit}}}} = {rg_prefactor:.2f} N^{{{rg_exponent:.2f}}}$')\n", + "plt.errorbar(N_MONOMERS, rg_summary[:,0],\n", + " yerr=rg_summary[:,1],\n", + " ls='', marker='o', capsize=5, capthick=1,\n", + " label=r'$R_g^{\\mathrm{simulation}}$')\n", + "plt.xlabel('Number of monomers $N$')\n", + "plt.ylabel('Radius of gyration [$\\sigma$]')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rf_summary[:,0]**2 / rg_summary[:,0]**2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the hydrodynamic radius $R_h$ of the polymers as a function of the number of monomers. What relation do you observe?" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The hydrodynamic radius can be calculated via the Stokes radius, i.e. the radius of a sphere that\n", + "diffuses at the same rate as the polymer. An approximative formula is $R_h \\approx c_h N^{1/3}$\n", + "with $c_h$ a constant." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "rh_summary = standard_error_mean_autocorrelation(rh_results, r'\\operatorname{acf}(R_h)')\n", + "rh_exponent, rh_prefactor = np.polyfit(np.log(N_MONOMERS), np.log(rh_summary[:,0]), 1)\n", + "rh_prefactor = np.exp(rh_prefactor)\n", + "\n", + "fig = plt.figure(figsize=(10, 6))\n", + "x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)\n", + "plt.plot(x, rh_prefactor * x**rh_exponent, '-',\n", + " label=rf'$R_h^{{\\mathrm{{fit}}}} = {rh_prefactor:.2f} N^{{{rh_exponent:.2f}}}$')\n", + "plt.errorbar(N_MONOMERS, rh_summary[:,0],\n", + " yerr=rh_summary[:,1],\n", + " ls='', marker='o', capsize=5, capthick=1,\n", + " label=r'$R_h^{\\mathrm{simulation}}$')\n", + "plt.xlabel('Number of monomers $N$')\n", + "plt.ylabel('Hydrodynamic radius [$\\sigma$]')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 3.2 Diffusion coefficient using the MSD method\n", + "\n", + "Calculate the diffusion coefficient of the polymers using the mean-squared displacement.\n", "\n", - "logging.info(\"Warming up the polymer chain.\")\n", - "system.time_step = 0.002\n", - "system.integrator.set_steepest_descent(f_max=1.0, gamma=10,\n", - " max_displacement=0.01)\n", - "system.integrator.run(2000)\n", - "system.integrator.set_vv()\n", - "logging.info(\"Warmup finished.\")\n", + "Recalling that for large $t$ the diffusion coefficient can be expressed as:\n", "\n", - "logging.info(\"Equilibration.\")\n", - "system.time_step = TIME_STEP\n", - "system.thermostat.set_langevin(kT=1.0, gamma=10, seed=42)\n", - "system.integrator.run(50000)\n", - "logging.info(\"Equilibration finished.\")\n", + "$$6D = \\lim_{t\\to\\infty} \\frac{\\partial \\operatorname{MSD}(t)}{\\partial t}$$\n", "\n", - "system.thermostat.turn_off()\n", + "which is simply the slope of the MSD in the diffusive mode." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# cutoff for the diffusive regime (approximative)\n", + "tau_f_index = 40\n", + "# cutoff for the data series (larger lag times have larger variance due to undersampling)\n", + "tau_max_index = 81\n", + "\n", + "plt.figure(figsize=(10, 10))\n", + "plt.xlabel(r'$\\tau$ [$\\Delta t$]')\n", + "plt.ylabel(r'MSD [$\\sigma^2$]')\n", + "for index, (tau, msd) in enumerate(zip(com_pos_tau_results, com_pos_msd_results)):\n", + " plt.loglog(tau[1:120], msd[1:120], label=f'N={N_MONOMERS[index]}')\n", + "plt.loglog(2 * [tau[tau_f_index]], [0, np.max(com_pos_msd_results)], '-', color='black')\n", + "plt.text(tau[tau_f_index], np.max(com_pos_msd_results), r'$\\tau_{f}$')\n", + "plt.loglog(2 * [tau[tau_max_index]], [0, np.max(com_pos_msd_results)], '-', color='black')\n", + "plt.text(tau[tau_max_index], np.max(com_pos_msd_results), r'$\\tau_{max}$')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "diffusion_msd = np.zeros(len(N_MONOMERS))\n", + "plt.figure(figsize=(10, 8))\n", + "plt.xlabel(r'$\\tau$ [$\\Delta t$]')\n", + "plt.ylabel(r'MSD [$\\sigma^2$]')\n", + "for index, (tau, msd) in enumerate(zip(com_pos_tau_results, com_pos_msd_results)):\n", + " a, b = np.polyfit(tau[tau_f_index:tau_max_index], msd[tau_f_index:tau_max_index], 1)\n", + " x = np.array([tau[1], tau[tau_max_index - 1]])\n", + " p = plt.plot(x, a * x + b, '-')\n", + " plt.plot(tau[1:tau_max_index], msd[1:tau_max_index], 'o', color=p[0].get_color(),\n", + " label=r'$N=${}'.format(N_MONOMERS[index]))\n", + " diffusion_msd[index] = a / 6\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the dependence of the diffusion coefficient on the hydrodynamic radius.\n", "\n", - "lbf = espressomd.lb.LBFluidGPU(kT=1, seed=123, agrid=1, dens=1, visc=5, tau=TIME_STEP)\n", - "system.actors.add(lbf)\n", - "system.thermostat.set_lb(LB_fluid=lbf, seed=142, gamma=5)\n", + "Recalling the formula for the diffusion coefficient of a short polymer in the Kirkwood–Zimm model:\n", "\n", - "logging.info(\"Warming up the system with LB fluid.\")\n", - "system.integrator.run(1000)\n", - "logging.info(\"LB fluid warming finished.\")\n", + "$$D = \\frac{D_0}{N} + \\frac{k_B T}{6 \\pi \\eta} \\left\\langle \\frac{1}{R_h} \\right\\rangle$$\n", "\n", - "# configure correlator\n", - "com_pos = espressomd.observables.ComPosition(ids=range(N_MONOMERS))\n", - "correlator = espressomd.accumulators.Correlator(obs1=com_pos, tau_lin=16, tau_max=LOOPS * STEPS, delta_N=1,\n", - " corr_operation=\"square_distance_componentwise\", compress1=\"discard1\")\n", - "system.auto_update_accumulators.add(correlator)\n", + "where $D_0$ is the monomer diffusion coefficient and $\\eta$ is the viscosity of the fluid.\n", + "For the Rouse model (implicit solvent), the second term disappears.\n", "\n", - "logging.info(\"Sampling started.\")\n", - "for i in range(LOOPS):\n", - " system.integrator.run(STEPS)\n", + "Hint: use $D = \\alpha N^{-1}$ and solve for $\\alpha$ in the Rouse model,\n", + "or use $D = \\alpha_1 N^{-1} + \\alpha_2 N^{-\\beta}$ with `rh_exponent` for $\\beta$\n", + "and solve for $\\alpha_1, \\alpha_2$ in the Zimm model." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import scipy.optimize\n", + "\n", + "def fitting_polymer_theory(polymer_model, n_monomers, diffusion, rh_exponent):\n", + " def rouse(x, a):\n", + " return a / x\n", + "\n", + " def kirkwood_zimm(x, a, b, exponent):\n", + " return a / x + b / x**exponent\n", + "\n", + " x = np.linspace(min(n_monomers) - 0.5, max(n_monomers) + 0.5, 20)\n", + "\n", + " if polymer_model == 'Rouse':\n", + " (a,), _ = scipy.optimize.curve_fit(rouse, n_monomers, diffusion)\n", + " label = rf'$D^{{\\mathrm{{fit}}}} = \\frac{{{a:.2f}}}{{N}}$'\n", + " y = rouse(x, a)\n", + " elif polymer_model == 'Zimm':\n", + " fitting_function = kirkwood_zimm\n", + " (a, b), _ = scipy.optimize.curve_fit(\n", + " lambda x, a, b: kirkwood_zimm(x, a, b, rh_exponent), n_monomers, diffusion)\n", + " y = kirkwood_zimm(x, a, b, rh_exponent)\n", + " label = f'''\\\n", + " $D^{{\\\\mathrm{{fit}}}} = \\\n", + " \\\\frac{{{a:.2f}}}{{N}} + \\\n", + " \\\\frac{{{b * 6 * np.pi:.3f} }}{{6\\\\pi}} \\\\cdot \\\n", + " \\\\frac{{{1}}}{{N^{{{rh_exponent:.2f}}}}}$ \\\n", + " '''\n", + " return x, y, label\n", + "\n", + "fig = plt.figure(figsize=(10, 6))\n", + "x, y, label = fitting_polymer_theory(POLYMER_MODEL, N_MONOMERS, diffusion_msd, rh_exponent)\n", + "plt.plot(x, y, '-', label=label)\n", + "plt.plot(N_MONOMERS, diffusion_msd, 'o', label=r'$D^{\\mathrm{simulation}}$')\n", + "plt.xlabel('Number of monomers $N$')\n", + "plt.ylabel(r'Diffusion coefficient [$\\sigma^2/t$]')\n", + "plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 0))\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 3.3 Diffusion coefficient using the Green–Kubo method\n", "\n", - "logging.info(\"Sampling finished.\")\n", + "Plot the autocorrelation function and check that the decay is roughly exponential.\n", "\n", - "correlator.finalize()\n", - "corrdata = correlator.result()\n", + "Hint: use $D = \\alpha e^{-\\beta \\tau}$ and solve for $\\alpha, \\beta$. You can leave out\n", + "the first data point in the ACF if necessary, and limit the fit to the stable region\n", + "in the first 20 data points." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def exponential(x, a, b):\n", + " return a * np.exp(-b * x)\n", + "\n", + "fig = plt.figure(figsize=(10, 8))\n", + "for N, tau, acf in zip(N_MONOMERS, com_vel_tau_results, com_vel_acf_results):\n", + " popt, _ = scipy.optimize.curve_fit(exponential, tau[:20], acf[:20])\n", + " x = np.linspace(tau[0], tau[20-1], 100)\n", + " p = plt.plot(x, exponential(x, *popt), '-')\n", + " plt.plot(tau[:20], acf[:20], 'o',\n", + " color=p[0].get_color(), label=rf'$R(\\tau)$ for N = {N}')\n", + "plt.xlabel(r'$\\tau$')\n", + "plt.ylabel('Autocorrelation function')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Green–Kubo integral for the diffusion coefficient take the following form:\n", "\n", - "tau = correlator.lag_times()\n", - "msd = corrdata[:, 0] + corrdata[:, 1] + corrdata[:, 2]\n", + "$$D = \\frac{1}{3} \\int_0^{+\\infty} \\left<\\vec{v_c}(\\tau)\\cdot\\vec{v_c}(0)\\right>\\, \\mathrm{d}\\tau$$\n", "\n", - "rh = system.analysis.calc_rh(chain_start=0, number_of_chains=1, chain_length=N_MONOMERS)[0]" + "Since our simulation is finite in time, we need to integrate up until $\\tau_{\\mathrm{int}}$. To find\n", + "the optimal value of $\\tau_{\\mathrm{int}}$, plot the integral as a function of $\\tau_{\\mathrm{int}}$\n", + "until you see a plateau. This plateau is usually followed by strong oscillations due to low\n", + "statistics in the long time tail of the autocorrelation function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "diffusion_gk = []\n", + "fig = plt.figure(figsize=(10, 6))\n", + "for N, tau, acf in zip(N_MONOMERS, com_vel_tau_results, com_vel_acf_results):\n", + " x = np.arange(2, 28)\n", + " y = [1 / 3 * np.trapz(acf[:j], tau[:j]) for j in x]\n", + " plt.plot(tau[x], y, label=rf'$D(\\tau_{{\\mathrm{{int}}}})$ for $N = {N}$')\n", + " diffusion_gk.append(np.mean(y[10:]))\n", + "plt.xlabel(r'$\\tau_{\\mathrm{int}}$')\n", + "plt.ylabel(r'$\\frac{1}{3} \\int_{\\tau=0}^{\\tau_{\\mathrm{int}}} \\left<\\vec{v_c}(\\tau)\\cdot\\vec{v_c}(0)\\right>\\, \\mathrm{d}\\tau$')\n", + "plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 0))\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot the dependence of the diffusion coefficient on the hydrodynamic radius." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig = plt.figure(figsize=(10, 8))\n", + "x, y, label = fitting_polymer_theory(POLYMER_MODEL, N_MONOMERS, diffusion_gk, rh_exponent)\n", + "plt.plot(x, y, '-', label=label)\n", + "plt.plot(N_MONOMERS, diffusion_gk, 'o', label=r'$D^{\\mathrm{simulation}}$')\n", + "plt.xlabel('Number of monomers $N$')\n", + "plt.ylabel(r'Diffusion coefficient [$\\sigma^2/t$]')\n", + "plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 0))\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us compare the value of the diffusion coefficients calcualted with the MSD and Green–Kubo methods:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(f'N\\tMSD\\t\\tGK\\t\\tdifference')\n", + "for N, d_msd, d_gk in zip(N_MONOMERS, diffusion_msd, diffusion_gk):\n", + " print(f'{N}\\t{d_msd:.2e}\\t{d_gk:.2e}\\t{np.ceil(np.abs(d_msd-d_gk) * 100 / d_msd):.0f}%')" ] }, { @@ -166,7 +791,9 @@ "[4] P. G. de Gennes. *Scaling Concepts in Polymer Physics*. Cornell University Press, Ithaca, NY, 1979. \n", "[5] M. Doi. *Introduction to Polymer Physics.* Clarendon Press, Oxford, 1996. \n", "[6] Michael Rubinstein and Ralph H. Colby. *Polymer Physics.* Oxford University Press, Oxford, UK, 2003. \n", - "[7] Daan Frenkel and Berend Smit. *Understanding Molecular Simulation.* Academic Press, San Diego, second edition, 2002." + "[7] Daan Frenkel and Berend Smit. *Understanding Molecular Simulation.* Academic Press, San Diego, second edition, 2002. \n", + "[8] W. Janke, Statistical analysis of simulations: Data correlations and error estimation, *Quantum Simulations of Complex Many-Body Systems: From Theory to Algorithms*, Lecture Notes, J. Grotendorst, D. Marx, A. Muramatsu (Eds.), John von Neumann Institute for Computing, 10:423–445, 2002. https://www.physik.uni-leipzig.de/~janke/Paper/nic10_423_2002.pdf \n", + "[9] M. Weigel, W. Janke, Error estimation and reduction with cross correlations, *Phys. Rev. E*, 81:066701, 2010, doi:[10.1103/PhysRevE.81.066701](https://doi.org/10.1103/PhysRevE.81.066701); Erratum-ibid 81:069902, 2010, doi:[10.1103/PhysRevE.81.069902](https://doi.org/10.1103/PhysRevE.81.069902). " ] } ], diff --git a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part4.ipynb b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part4.ipynb index 460e79ac282..6ec60707b67 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part4.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part4.ipynb @@ -11,7 +11,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### 6 Poiseuille flow in ESPResSo\n", + "## Poiseuille flow in ESPResSo\n", "\n", "Poiseuille flow is the flow through a pipe or (in our case) a slit\n", "under a homogeneous force density, e.g. gravity. In the limit of small Reynolds\n", @@ -41,6 +41,13 @@ "applied to every node." ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1. Setting up the system" + ] + }, { "cell_type": "code", "execution_count": null, @@ -64,41 +71,72 @@ "\n", "# System constants\n", "BOX_L = 16.0\n", + "TIME_STEP = 0.01\n", + "\n", + "system = espressomd.System(box_l=[BOX_L]*3)\n", + "system.time_step = TIME_STEP\n", + "system.cell_system.skin = 0.4" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 1.1 Setting up the lattice-Boltzmann fluid\n", + "\n", + "We will now create a lattice-Boltzmann fluid confined between two walls." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# LB parameters\n", "AGRID = 0.5\n", "VISCOSITY = 2.0\n", "FORCE_DENSITY = [0.0, 0.001, 0.0]\n", "DENSITY = 1.5\n", - "TIME_STEP = 0.01\n", "\n", - "system = espressomd.System(box_l=[BOX_L]*3)\n", - "system.time_step = TIME_STEP\n", - "system.cell_system.skin = 0.4\n", + "# LB boundary parameters\n", + "WALL_OFFSET = AGRID" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden", + "solution2_first": true + }, + "source": [ + "Create a lattice-Boltzmann actor and append it to the list of system actors. Use the GPU implementation of LB.\n", "\n", + "You can refer to section [setting up a LB fluid](http://espressomd.org/html/doc/lb.html#setting-up-a-lb-fluid)\n", + "in the user guide." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "solution2": "hidden" + }, + "source": [ + "```python\n", + "logging.info(\"Setup LB fluid.\")\n", "lbf = espressomd.lb.LBFluidGPU(agrid=AGRID, dens=DENSITY, visc=VISCOSITY, tau=TIME_STEP,\n", - " ext_force_density=FORCE_DENSITY)\n", + " ext_force_density=FORCE_DENSITY)\n", "system.actors.add(lbf)\n", - "\n", - "logging.info(\"Setup LB boundaries.\")\n", - "WALL_OFFSET = AGRID\n", - "top_wall = espressomd.lbboundaries.LBBoundary(shape=espressomd.shapes.Wall(normal=[1, 0, 0], dist=WALL_OFFSET))\n", - "bottom_wall = espressomd.lbboundaries.LBBoundary(shape=espressomd.shapes.Wall(normal=[-1, 0, 0], dist=-(BOX_L - WALL_OFFSET)))\n", - "\n", - "system.lbboundaries.add(top_wall)\n", - "system.lbboundaries.add(bottom_wall)\n", - "\n", - "logging.info(\"Iterate until the flow profile converges (5000 LB updates).\")\n", - "system.integrator.run(5000)\n", - "\n", - "logging.info(\"Extract fluid velocity along the x-axis\")\n", - "fluid_velocities = np.zeros(lbf.shape[0])\n", - "for x in range(lbf.shape[0]):\n", - " # Average over the nodes in y direction\n", - " v_tmp = np.zeros(lbf.shape[1])\n", - " for y in range(lbf.shape[1]):\n", - " v_tmp[y] = lbf[x, y, 0].velocity[1]\n", - " fluid_velocities[x] = np.average(v_tmp)" + "```" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "metadata": { @@ -106,7 +144,9 @@ "solution2_first": true }, "source": [ - "Use the data to fit a parabolic function. Can you confirm the analytic solution?" + "Create a LB boundary and append it to the list of system LB boundaries.\n", + "\n", + "You can refer to section [using shapes as lattice-Boltzmann boundary](http://espressomd.org/html/doc/lb.html#using-shapes-as-lattice-boltzmann-boundary) in the user guide." ] }, { @@ -116,39 +156,89 @@ }, "source": [ "```python\n", - "# Extract fluid velocity along the x-axis\n", - "fluid_velocities = np.zeros((lbf.shape[0], 2))\n", + "logging.info(\"Setup LB boundaries.\")\n", + "top_wall = espressomd.shapes.Wall(normal=[1, 0, 0], dist=WALL_OFFSET)\n", + "bottom_wall = espressomd.shapes.Wall(normal=[-1, 0, 0], dist=-(BOX_L - WALL_OFFSET))\n", + "\n", + "top_boundary = espressomd.lbboundaries.LBBoundary(shape=top_wall)\n", + "bottom_boundary = espressomd.lbboundaries.LBBoundary(shape=bottom_wall)\n", + "\n", + "system.lbboundaries.add(top_boundary)\n", + "system.lbboundaries.add(bottom_boundary)\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2. Simulation\n", + "\n", + "We will now simulate the fluid flow until we reach the steady state." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "logging.info(\"Iterate until the flow profile converges (5000 LB updates).\")\n", + "system.integrator.run(5000)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 3. Data analysis\n", + "\n", + "We can now extract the flow profile and compare it to the analytical solution for the planar Poiseuille flow." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "logging.info(\"Extract fluid velocities along the x-axis\")\n", + "fluid_positions = np.zeros(lbf.shape[0])\n", + "fluid_velocities = np.zeros(lbf.shape[0])\n", "for x in range(lbf.shape[0]):\n", - " # Average over the node in y direction\n", - " v_tmp = np.zeros(lbf.shape[1])\n", - " for y in range(lbf.shape[1]):\n", - " v_tmp[y] = lbf[x, y, 0].velocity[1]\n", - " fluid_velocities[x, 0] = (x + 0.5) * AGRID\n", - " fluid_velocities[x, 1] = np.average(v_tmp)\n", + " # Average over the nodes in y direction\n", + " v = [lbf[x, y, 0].velocity[1] for y in range(lbf.shape[1])]\n", + " fluid_velocities[x] = np.average(v)\n", + " fluid_positions[x] = (x + 0.5) * AGRID\n", "\n", "def poiseuille_flow(x, force_density, dynamic_viscosity, height):\n", - " return force_density / (2.0 * dynamic_viscosity) * \\\n", - " (height**2.0 / 4.0 - x**2.0)\n", + " return force_density / (2 * dynamic_viscosity) * (height**2 / 4 - x**2)\n", "\n", "# Note that the LB viscosity is not the dynamic viscosity but the\n", "# kinematic viscosity (mu=LB_viscosity * density)\n", "x_values = np.linspace(0.0, BOX_L, lbf.shape[0])\n", "HEIGHT = BOX_L - 2.0 * AGRID\n", "# analytical curve\n", - "y_values = poiseuille_flow(x_values - (HEIGHT / 2.0 + AGRID), FORCE_DENSITY[1],\n", + "y_values = poiseuille_flow(x_values - (HEIGHT / 2 + AGRID), FORCE_DENSITY[1],\n", " VISCOSITY * DENSITY, HEIGHT)\n", - "# velocity is zero outside the walls\n", + "# velocity is zero inside the walls\n", "y_values[np.nonzero(x_values < WALL_OFFSET)] = 0.0\n", "y_values[np.nonzero(x_values > BOX_L - WALL_OFFSET)] = 0.0\n", "\n", - "fig1 = plt.figure(num=None, figsize=(10, 6), dpi=80, facecolor='w', edgecolor='k')\n", + "fig1 = plt.figure(figsize=(10, 6))\n", "plt.plot(x_values, y_values, '-', linewidth=2, label='analytical')\n", - "plt.plot(fluid_velocities[:, 0], fluid_velocities[:, 1], 'o', label='simulation')\n", + "plt.plot(fluid_positions, fluid_velocities, 'o', label='simulation')\n", "plt.xlabel('Position on the $x$-axis', fontsize=16)\n", "plt.ylabel('Fluid velocity in $y$-direction', fontsize=16)\n", "plt.legend()\n", - "plt.show()\n", - "```" + "plt.show()" ] }, { diff --git a/doc/tutorials/04-lattice_boltzmann/CMakeLists.txt b/doc/tutorials/04-lattice_boltzmann/CMakeLists.txt index 85f800c7f72..7777020c366 100644 --- a/doc/tutorials/04-lattice_boltzmann/CMakeLists.txt +++ b/doc/tutorials/04-lattice_boltzmann/CMakeLists.txt @@ -2,8 +2,7 @@ configure_tutorial_target( TARGET tutorial_04 DEPENDS 04-lattice_boltzmann_part1.ipynb 04-lattice_boltzmann_part2.ipynb 04-lattice_boltzmann_part3.ipynb 04-lattice_boltzmann_part4.ipynb figures/latticeboltzmann-grid.png - figures/latticeboltzmann-momentumexchange.png - scripts/04-lattice_boltzmann_part3_solution.py) + figures/latticeboltzmann-momentumexchange.png) nb_export(TARGET tutorial_04 SUFFIX "1" FILE "04-lattice_boltzmann_part1.ipynb" HTML_RUN) diff --git a/doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part3_solution.py b/doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part3_solution.py deleted file mode 100644 index cfcc927465d..00000000000 --- a/doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part3_solution.py +++ /dev/null @@ -1,125 +0,0 @@ -# Copyright (C) 2019 The ESPResSo project -# -# This file is part of ESPResSo. -# -# ESPResSo is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# ESPResSo is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - -import logging -import numpy as np -logging.basicConfig(level=logging.INFO) - -import espressomd -espressomd.assert_features(['LENNARD_JONES']) -import espressomd.accumulators -import espressomd.observables -import espressomd.polymer - -# Setup constant -TIME_STEP = 0.01 -LOOPS = 5000 -STEPS = 100 - -# System setup -system = espressomd.System(box_l=[32.0, 32.0, 32.0]) -system.cell_system.skin = 0.4 - -# Lennard-Jones interaction -system.non_bonded_inter[0, 0].lennard_jones.set_params( - epsilon=1.0, sigma=1.0, - shift="auto", cutoff=2.0**(1.0 / 6.0)) - -# Fene interaction -fene = espressomd.interactions.FeneBond(k=7, d_r_max=2) -system.bonded_inter.add(fene) - -N_MONOMERS = [10, 20, 40] - -rh_results = np.zeros(len(N_MONOMERS)) -diffusion_constant_results = np.zeros(len(N_MONOMERS)) -for index, N in enumerate(N_MONOMERS): - logging.info("Polymer size: {}".format(N)) - system.part.clear() - system.thermostat.turn_off() - system.actors.clear() - system.auto_update_accumulators.clear() - # Setup polymer of part_id 0 with fene bond - positions = espressomd.polymer.linear_polymer_positions(n_polymers=1, - beads_per_chain=N, - bond_length=1, seed=5642, - min_distance=0.9) - for i, pos in enumerate(positions[0]): - pid = len(system.part) - system.part.add(id=pid, pos=pos) - if i > 0: - system.part[pid].add_bond((fene, pid - 1)) - - logging.info("Warming up the polymer chain.") - system.time_step = 0.002 - system.integrator.set_steepest_descent( - f_max=1.0, - gamma=10, - max_displacement=0.01) - system.integrator.run(2000) - system.integrator.set_vv() - logging.info("Warmup finished.") - - logging.info("Equilibration.") - system.time_step = TIME_STEP - system.thermostat.set_langevin(kT=1.0, gamma=10, seed=42) - system.integrator.run(50000) - logging.info("Equilibration finished.") - - system.thermostat.turn_off() - - lbf = espressomd.lb.LBFluidGPU( - kT=1, - seed=123, - agrid=1, - dens=1, - visc=5, - tau=TIME_STEP) - system.actors.add(lbf) - system.thermostat.set_lb(LB_fluid=lbf, seed=142, gamma=5) - - logging.info("Warming up the system with LB fluid.") - system.integrator.run(1000) - logging.info("Warming up the system with LB fluid finished.") - - # configure correlator - com_pos = espressomd.observables.ComPosition(ids=range(N)) - correlator = espressomd.accumulators.Correlator( - obs1=com_pos, tau_lin=16, tau_max=LOOPS * STEPS, delta_N=1, - corr_operation="square_distance_componentwise", compress1="discard1") - system.auto_update_accumulators.add(correlator) - - logging.info("Sampling started.") - rhs = np.zeros(LOOPS) - for i in range(LOOPS): - system.integrator.run(STEPS) - rhs[i] = system.analysis.calc_rh( - chain_start=0, - number_of_chains=1, - chain_length=N)[0] - - logging.info("Sampling finished.") - - correlator.finalize() - corrdata = correlator.result() - - rh_results[index] = np.average(rhs) - tau = correlator.lag_times()[1:] - msd = corrdata[1:, 0] + corrdata[1:, 1] + corrdata[1:, 2] - np.save('msd_{}'.format(N), np.c_[tau, msd]) - -np.save('rh.npy', rh_results) diff --git a/doc/tutorials/04-lattice_boltzmann/tutor_notes.md b/doc/tutorials/04-lattice_boltzmann/tutor_notes.md new file mode 100644 index 00000000000..ac7a0576080 --- /dev/null +++ b/doc/tutorials/04-lattice_boltzmann/tutor_notes.md @@ -0,0 +1,45 @@ +# Part 1: Lattice-Boltzmann + +## Physical learning goals + +* Learning basic concepts of the LBM + +## Espresso learning goals + +* Learning the LB interface in ESPResSo + +# Part 2: Brownian motion of a sphere in a fluid + +## Physical learning goals + +* Getting a short introduction to the Flory theory of polymers +* Learning the basics of Brownian motion +* Modeling the Brownian motion of a spherical object +* Calculating the diffusion coefficient of a spherical object + +# Part 3: Brownian motion of a polymer in a fluid + +## Physical learning goals + +* Modeling the Brownian motion of a polymer +* Learning about the hydrodynamic radius, radius of gyration, end-to-end distance +* Reproducing theoretical results of the Flory theory of polymers +* Using the Green-Kubo relation to calculate the diffusion coefficient of a polymer +* Estimating the correlation-corrected standard error of the mean of a time series +* Observing the effect of implicit hydrodynamics vs. LB hydrodynamics on polymer diffusion + +## Espresso learning goals + +* Creating linear polymers with ESPResSo functions +* Using tools from the analysis module not based on the Observable framework +* Setting up correlators for various operations, with compression + +# Part 4: Planar Poiseuille flow + +## Physical learning goals + +* Modeling a laminar fluid flow between two walls + +## Espresso learning goals + +* Setting up an LB fluid and LB boundaries diff --git a/src/python/espressomd/lb.pyx b/src/python/espressomd/lb.pyx index 9fd7697e50d..e1bd36ba0e0 100644 --- a/src/python/espressomd/lb.pyx +++ b/src/python/espressomd/lb.pyx @@ -44,6 +44,36 @@ def assert_agrid_tau_set(obj): cdef class HydrodynamicInteraction(Actor): + """ + Base class for LB implementations. + + Parameters + ---------- + agrid : :obj:`float` + Lattice constant. The box size in every direction must be an integer + multiple of ``agrid``. + tau : :obj:`float` + LB time step. The MD time step must be an integer multiple of ``tau``. + dens : :obj:`float` + Fluid density. + visc : :obj:`float` + Fluid kinematic viscosity. + bulk_visc : :obj:`float`, optional + Fluid bulk viscosity. + gamma_odd : :obj:`int`, optional + Relaxation parameter :math:`\\gamma_{\\textrm{odd}}` for kinetic modes. + gamma_even : :obj:`int`, optional + Relaxation parameter :math:`\\gamma_{\\textrm{even}}` for kinetic modes. + ext_force_density : (3,) array_like of :obj:`float`, optional + Force density applied on the fluid. + kT : :obj:`float`, optional + Thermal energy of the simulated heat bath (for thermalized fluids). + Set it to 0 for an unthermalized fluid. + seed : :obj:`int`, optional + Initial counter value (or seed) of the philox RNG. + Required for a thermalized fluid. Must be positive. + """ + def _lb_init(self): raise Exception( "Subclasses of HydrodynamicInteraction must define the _lb_init() method.") @@ -143,13 +173,15 @@ cdef class HydrodynamicInteraction(Actor): Parameters ---------- - interpolation_order : :obj:`str` - ``"linear"`` for linear interpolation, ``"quadratic"`` for quadratic interpolation. + interpolation_order : :obj:`str`, \{"linear", "quadratic"\} + ``"linear"`` for trilinear interpolation, ``"quadratic"`` for + quadratic interpolation. For the CPU implementation of LB, only + ``"linear"`` is available. """ - if (interpolation_order == "linear"): + if interpolation_order == "linear": lb_lbinterpolation_set_interpolation_order(linear) - elif (interpolation_order == "quadratic"): + elif interpolation_order == "quadratic": lb_lbinterpolation_set_interpolation_order(quadratic) else: raise ValueError("Invalid parameter") @@ -160,7 +192,7 @@ cdef class HydrodynamicInteraction(Actor): Parameters ---------- pos : (3,) array_like of :obj:`float` - The position at which velocity is requested. + The position at which velocity is requested. Returns ------- @@ -309,6 +341,7 @@ cdef class HydrodynamicInteraction(Actor): cdef class LBFluid(HydrodynamicInteraction): """ Initialize the lattice-Boltzmann method for hydrodynamic flow using the CPU. + See :class:`HydrodynamicInteraction` for the list of parameters. """ @@ -324,6 +357,7 @@ IF CUDA: cdef class LBFluidGPU(HydrodynamicInteraction): """ Initialize the lattice-Boltzmann method for hydrodynamic flow using the GPU. + See :class:`HydrodynamicInteraction` for the list of parameters. """ diff --git a/testsuite/scripts/tutorials/CMakeLists.txt b/testsuite/scripts/tutorials/CMakeLists.txt index 2a869227081..713013202ab 100644 --- a/testsuite/scripts/tutorials/CMakeLists.txt +++ b/testsuite/scripts/tutorials/CMakeLists.txt @@ -33,9 +33,9 @@ tutorial_test(FILE test_02-charged_system__scripts__nacl_units_confined.py) tutorial_test(FILE test_02-charged_system__scripts__nacl_units_confined_vis.py) tutorial_test(FILE test_02-charged_system__scripts__nacl_units.py) tutorial_test(FILE test_02-charged_system__scripts__nacl_units_vis.py) -tutorial_test(FILE test_04-lattice_boltzmann_part2.py LABELS "gpu") -tutorial_test(FILE test_04-lattice_boltzmann_part3.py LABELS "gpu") -tutorial_test(FILE test_04-lattice_boltzmann_part3_solution.py LABELS "gpu") +tutorial_test(FILE test_04-lattice_boltzmann_part2.py) +tutorial_test(FILE test_04-lattice_boltzmann_part3.py SUFFIX rouse) +tutorial_test(FILE test_04-lattice_boltzmann_part3.py SUFFIX zimm LABELS "gpu") tutorial_test(FILE test_04-lattice_boltzmann_part4.py LABELS "gpu") tutorial_test(FILE test_05-raspberry_electrophoresis.py LABELS "gpu") tutorial_test(FILE test_06-active_matter__flow_field.py LABELS "gpu") diff --git a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py index cee46943cc0..805b5445459 100644 --- a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py +++ b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py @@ -17,16 +17,31 @@ import unittest as ut import importlib_wrapper +import numpy as np +import scipy.optimize tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( - "@TUTORIALS_DIR@/04-lattice_boltzmann/04-lattice_boltzmann_part2.py", - gpu=True, LOOPS=400) + "@TUTORIALS_DIR@/04-lattice_boltzmann/04-lattice_boltzmann_part2.py") @skipIfMissingFeatures class Tutorial(ut.TestCase): system = tutorial.system + def test_ballistic_regime(self): + for (tau_p, tau, msd) in zip(tutorial.tau_p_values, + tutorial.tau_results, + tutorial.msd_results): + popt, _ = scipy.optimize.curve_fit( + tutorial.quadratic, tau[:tau_p], msd[:tau_p]) + residuals = msd[:tau_p] - tutorial.quadratic(tau[:tau_p], *popt) + np.testing.assert_allclose(residuals, 0, rtol=0, atol=1e-3) + + def test_diffusion_coefficient(self): + D_val = tutorial.diffusion_results + D_ref = tutorial.KT / np.array(tutorial.gammas) + np.testing.assert_allclose(D_val, D_ref, rtol=0, atol=0.1) + if __name__ == "__main__": ut.main() diff --git a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py index f7ad4624584..56075de1c3e 100644 --- a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py +++ b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py @@ -17,16 +17,39 @@ import unittest as ut import importlib_wrapper +import numpy as np + +if '@TEST_SUFFIX@' == 'rouse': + params = {} +elif '@TEST_SUFFIX@' == 'zimm': + params = {'LOOPS': 2000, 'POLYMER_MODEL': 'Zimm', 'gpu': True} tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( "@TUTORIALS_DIR@/04-lattice_boltzmann/04-lattice_boltzmann_part3.py", - gpu=True) + script_suffix="@TEST_SUFFIX@", **params) @skipIfMissingFeatures class Tutorial(ut.TestCase): system = tutorial.system + def test_exponents(self): + msg = 'The R_F exponent should be close to 0.588' + self.assertGreater(tutorial.rf_exponent, 0.50, msg=msg) + self.assertLess(tutorial.rf_exponent, 0.85, msg=msg) + msg = 'The R_g exponent should be close to 0.588' + self.assertGreater(tutorial.rg_exponent, 0.50, msg=msg) + self.assertLess(tutorial.rg_exponent, 0.75, msg=msg) + msg = 'The R_h exponent should be close to 0.333' + self.assertGreater(tutorial.rh_exponent, 0.30, msg=msg) + self.assertLess(tutorial.rh_exponent, 0.50, msg=msg) + + def test_diffusion_coefficients(self): + reference = [0.0363, 0.0269, 0.0234] + np.testing.assert_allclose( + tutorial.diffusion_msd, reference, rtol=0.15) + np.testing.assert_allclose(tutorial.diffusion_gk, reference, rtol=0.15) + if __name__ == "__main__": ut.main() diff --git a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3_solution.py b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3_solution.py deleted file mode 100644 index 81bcb0b7093..00000000000 --- a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3_solution.py +++ /dev/null @@ -1,33 +0,0 @@ -# Copyright (C) 2019 The ESPResSo project -# -# This file is part of ESPResSo. -# -# ESPResSo is free software: you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation, either version 3 of the License, or -# (at your option) any later version. -# -# ESPResSo is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - -import unittest as ut -import importlib_wrapper - - -tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( - "@TUTORIALS_DIR@/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part3_solution.py", - gpu=True, LOOPS=100, STEPS=1, N_MONOMERS=[10]) - - -@skipIfMissingFeatures -class Tutorial(ut.TestCase): - system = tutorial.system - - -if __name__ == "__main__": - ut.main() diff --git a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part4.py b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part4.py index 58dded8160c..eb2706efb8d 100644 --- a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part4.py +++ b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part4.py @@ -31,7 +31,7 @@ class Tutorial(ut.TestCase): def test_flow_profile(self): analytical = tutorial.y_values - simulation = tutorial.fluid_velocities[:, 1] + simulation = tutorial.fluid_velocities rmsd = np.sqrt(np.mean(np.square(analytical - simulation))) self.assertLess(rmsd, 2e-5 * tutorial.AGRID / tutorial.lbf.tau)