From 808e94fd0cbee51e2d841228f03d7d4f8c892e07 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 11 Sep 2020 16:53:29 +0200 Subject: [PATCH 01/16] Tutorial 04: Spell check and grammar check --- .../04-lattice_boltzmann_part1.ipynb | 11 ++++++----- .../04-lattice_boltzmann_part2.ipynb | 8 ++++---- 2 files changed, 10 insertions(+), 9 deletions(-) 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..c82aa0ecfa3 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb @@ -13,7 +13,8 @@ "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." @@ -96,7 +97,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", @@ -197,7 +198,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", @@ -314,14 +315,14 @@ "### 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", + "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 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. One way to construct a wall is:\n", "\n", "```python\n", "import espressomd.lbboundaries\n", 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..e06bf238e5c 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb @@ -29,7 +29,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 +37,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 +51,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", @@ -118,7 +118,7 @@ "# Setup observable\n", "pos = espressomd.observables.ParticlePositions(ids=(0,))\n", "\n", - "# Run for different friction constants\n", + "# Run for different friction coefficients\n", "lb_gammas = [1.0, 2.0, 4.0, 10.0]\n", "msd_results = []\n", "\n", From bed58d97eedd1ae2117e54154c9a62dbb14e82b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 11 Sep 2020 17:53:31 +0200 Subject: [PATCH 02/16] python: Document LB parameters --- .../04-lattice_boltzmann_part1.ipynb | 6 ++- src/python/espressomd/lb.pyx | 41 +++++++++++++++++-- 2 files changed, 41 insertions(+), 6 deletions(-) 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 c82aa0ecfa3..f64eaa2c344 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb @@ -200,7 +200,7 @@ "\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", + "via linear interpolation and a force is applied on the particle\n", "that is proportional to the velocity difference between particle and fluid:\n", "\n", "\\begin{equation}\n", @@ -214,7 +214,7 @@ "
\n", "\n", "
\n", - "
The coupling scheme between fluid and particles is based on the interpolation of the fluid velocity $\\vec{u}$ from the grid nodes. This is done by trilinear interpolation. The difference between the particle velocity $\\vec{v}(t)$ and the interpolated velocity $\\vec{u}(\\vec{r},t)$ is used in the momentum exchange of the equation $\\vec{F}_H$ above.
\n", + "
The coupling scheme between fluid and particles is based on the interpolation of the fluid velocity $\\vec{u}$ from the grid nodes. This is done by linear interpolation. The difference between the particle velocity $\\vec{v}(t)$ and the interpolated velocity $\\vec{u}(\\vec{r},t)$ is used in the momentum exchange of the equation $\\vec{F}_H$ above.
\n", "
\n", "
" ] @@ -269,6 +269,8 @@ "+ 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", diff --git a/src/python/espressomd/lb.pyx b/src/python/espressomd/lb.pyx index 9fd7697e50d..516a3d0c204 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,14 @@ cdef class HydrodynamicInteraction(Actor): Parameters ---------- - interpolation_order : :obj:`str` - ``"linear"`` for linear interpolation, ``"quadratic"`` for quadratic interpolation. + interpolation_order : :obj:`str`, \{"linear", "quadratic"\} + Interpolation order. 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") @@ -309,6 +340,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 +356,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. """ From c1aed230ddee76c6d5feee575af6f82b47b6987b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 11 Sep 2020 17:59:36 +0200 Subject: [PATCH 03/16] tutorial 04: part 1: Update text --- .../04-lattice_boltzmann_part1.ipynb | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) 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 f64eaa2c344..e810e43d428 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb @@ -9,7 +9,7 @@ "#### 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", @@ -47,7 +47,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", @@ -56,6 +56,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", @@ -150,7 +151,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", @@ -172,9 +173,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", @@ -321,10 +320,7 @@ "\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 instantiated by passing a shape object to the LBBoundary class. One way to construct 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", @@ -335,7 +331,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." ] From 4dfbd96091b176f5c05f19fbc780fcbf484387b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 11 Sep 2020 19:03:11 +0200 Subject: [PATCH 04/16] tutorial 04: part 2: Simplify code Store the lag time and MSD series separately. Remove the superfluous LOOPS variable. Indicate which section in Frenkel 2002 discusses the diffusion formula. --- .../04-lattice_boltzmann_part2.ipynb | 28 ++++++++----------- .../test_04-lattice_boltzmann_part2.py | 2 +- 2 files changed, 12 insertions(+), 18 deletions(-) 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 e06bf238e5c..345d5897e2e 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb @@ -98,21 +98,17 @@ "\n", "logging.basicConfig(level=logging.INFO, stream=sys.stdout)\n", "\n", - "\n", "# Constants\n", - "LOOPS = 40000\n", - "STEPS = 10\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", @@ -120,6 +116,7 @@ "\n", "# Run for different friction coefficients\n", "lb_gammas = [1.0, 2.0, 4.0, 10.0]\n", + "tau_results = []\n", "msd_results = []\n", "\n", "for gamma in lb_gammas:\n", @@ -131,17 +128,16 @@ " system.integrator.run(1000)\n", " logging.info(\"Equilibration finished.\")\n", " # Setup observable correlator\n", - " correlator = espressomd.accumulators.Correlator(obs1=pos, tau_lin=16, tau_max=LOOPS * STEPS, delta_N=1,\n", + " correlator = espressomd.accumulators.Correlator(\n", + " obs1=pos, tau_lin=16, tau_max=STEPS, delta_N=1,\n", " corr_operation=\"square_distance_componentwise\", compress1=\"discard1\")\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.\")" ] @@ -167,13 +163,11 @@ "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", + "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(lb_gammas[index]))\n", "plt.legend()\n", "plt.show()" ] @@ -202,7 +196,7 @@ "[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." ] } ], diff --git a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py index cee46943cc0..59e73d9927b 100644 --- a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py +++ b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py @@ -20,7 +20,7 @@ tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( "@TUTORIALS_DIR@/04-lattice_boltzmann/04-lattice_boltzmann_part2.py", - gpu=True, LOOPS=400) + gpu=True, STEPS=4000) @skipIfMissingFeatures From a9a94cf19ef82f97b1c63530fb1a019b899431cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 11 Sep 2020 22:56:25 +0200 Subject: [PATCH 05/16] tutorial 04: part 2: Provide solutions Add solutions and literature references. Test the solutions in CI. --- .../04-lattice_boltzmann_part2.ipynb | 143 +++++++++++++++++- .../test_04-lattice_boltzmann_part2.py | 18 ++- 2 files changed, 155 insertions(+), 6 deletions(-) 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 345d5897e2e..054918498e6 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb @@ -160,9 +160,9 @@ "\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", + "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", @@ -181,7 +181,138 @@ "random forces on the particle and within the LB fluid 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(lb_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(lb_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')\n", + "gammas = np.linspace(0.9 * min(lb_gammas), 1.1 * max(lb_gammas), 50)\n", + "plt.plot(gammas, lbf.kT / gammas, '-', label=r'$k_BT\\gamma^{-1}$')\n", + "plt.plot(lb_gammas, diffusion_results, 'o', label='D')\n", + "plt.legend()\n", + "plt.show()" ] }, { @@ -196,7 +327,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*, section 4.4.1. 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/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py index 59e73d9927b..d71cc5fa847 100644 --- a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py +++ b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py @@ -17,16 +17,32 @@ 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, STEPS=4000) + gpu=True, STEPS=8000, cutoff_limit=66) @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.lbf.kT / np.array(tutorial.lb_gammas) + np.testing.assert_allclose(D_val, D_ref, rtol=0, atol=0.1) + if __name__ == "__main__": ut.main() From 32c410f33bcf0e7bc9230d76e0ea5bc967de8eb8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Sat, 19 Sep 2020 22:59:34 +0200 Subject: [PATCH 06/16] docs: Fix documentation --- .../04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb | 4 ++-- .../04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb | 6 +++--- src/python/espressomd/lb.pyx | 5 +++-- 3 files changed, 8 insertions(+), 7 deletions(-) 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 e810e43d428..b9c7b6bc2b2 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb @@ -199,7 +199,7 @@ "\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 linear interpolation and a force is applied on the particle\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", "\n", "\\begin{equation}\n", @@ -213,7 +213,7 @@ "
\n", "\n", "
\n", - "
The coupling scheme between fluid and particles is based on the interpolation of the fluid velocity $\\vec{u}$ from the grid nodes. This is done by linear interpolation. The difference between the particle velocity $\\vec{v}(t)$ and the interpolated velocity $\\vec{u}(\\vec{r},t)$ is used in the momentum exchange of the equation $\\vec{F}_H$ above.
\n", + "
The coupling scheme between fluid and particles is based on the interpolation of the fluid velocity $\\vec{u}$ from the grid nodes. This is done by trilinear interpolation. The difference between the particle velocity $\\vec{v}(t)$ and the interpolated velocity $\\vec{u}(\\vec{r},t)$ is used in the momentum exchange of the equation $\\vec{F}_H$ above.
\n", "
\n", "
" ] 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 9e1b0d12940..798cd6dbe50 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb @@ -18,11 +18,11 @@ "```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", + " 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", + "$r_\\text{cut} = 1.1225$ and $\\varepsilon_\\text{shift}=0.25$ between particles\n", "of type 0, which is the desired \n", "repulsive interaction. The command\n", "\n", @@ -76,7 +76,7 @@ "\n", "espressomd.assert_features(['LENNARD_JONES'])\n", "\n", - "# Setup constant\n", + "# Setup constants\n", "TIME_STEP = 0.01\n", "LOOPS = 100\n", "STEPS = 100\n", diff --git a/src/python/espressomd/lb.pyx b/src/python/espressomd/lb.pyx index 516a3d0c204..e1bd36ba0e0 100644 --- a/src/python/espressomd/lb.pyx +++ b/src/python/espressomd/lb.pyx @@ -174,7 +174,8 @@ cdef class HydrodynamicInteraction(Actor): Parameters ---------- interpolation_order : :obj:`str`, \{"linear", "quadratic"\} - Interpolation order. For the CPU implementation of LB, only + ``"linear"`` for trilinear interpolation, ``"quadratic"`` for + quadratic interpolation. For the CPU implementation of LB, only ``"linear"`` is available. """ @@ -191,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 ------- From 216a2feef443f7d454b7d118fc078eec218a9d9c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Sat, 19 Sep 2020 23:06:20 +0200 Subject: [PATCH 07/16] tutorial 04: part 3: Move solution to the notebook Rewrite the solution to run faster with shorter polymers. Calculate several polymer statistics (R_F, R_g, R_h) and explain they depend on the number of monomers. --- .../04-lattice_boltzmann_part3.ipynb | 188 +++++++++++------- .../04-lattice_boltzmann/CMakeLists.txt | 3 +- .../04-lattice_boltzmann_part3_solution.py | 126 ------------ testsuite/scripts/tutorials/CMakeLists.txt | 1 - .../test_04-lattice_boltzmann_part3.py | 2 +- ...est_04-lattice_boltzmann_part3_solution.py | 33 --- 6 files changed, 119 insertions(+), 234 deletions(-) delete mode 100644 doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part3_solution.py delete mode 100644 testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3_solution.py 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 798cd6dbe50..deadc78d93b 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb @@ -17,17 +17,15 @@ "\n", "```python\n", "system.non_bonded_inter[0, 0].lennard_jones.set_params(\n", - " epsilon=1.0, sigma=1.0,\n", - " shift=0.25, cutoff=1.1225)\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.1225$ and $\\varepsilon_\\text{shift}=0.25$ between particles\n", - "of type 0, which is the desired \n", - "repulsive interaction. The command\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,20 @@ "```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", "\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", + "Furthermore we want to compute the diffusion constant of the polymer for different numbers of monomers.\n", + "For this purpose we can again use the multiple tau correlator. The following script computes the mean\n", + "squared displacement for the center of mass of the polymer as well as the average hydrodynamic radius\n", + "$R_h$, end-to-end distance $R_F$ and radius of gyration $R_g$.\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." + "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$." ] }, { @@ -70,7 +73,6 @@ "import espressomd.accumulators\n", "import espressomd.observables\n", "import espressomd.polymer\n", - "import espressomd.minimize_energy\n", "\n", "logging.basicConfig(level=logging.INFO, stream=sys.stdout)\n", "\n", @@ -78,79 +80,123 @@ "\n", "# Setup constants\n", "TIME_STEP = 0.01\n", - "LOOPS = 100\n", + "LOOPS = 4000\n", "STEPS = 100\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", - "# 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", - "\n", - "\n", - "logging.info(\"Warming up the polymer chain.\")\n", - "system.time_step = 0.002\n", - "espressomd.minimize_energy.steepest_descent(\n", - " system, f_max=1.0, gamma=10, max_steps=2000, max_displacement=0.01)\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=10, seed=42)\n", - "system.integrator.run(50000)\n", - "logging.info(\"Equilibration finished.\")\n", - "\n", - "system.thermostat.turn_off()\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", - "\n", - "logging.info(\"Warming up the system with LB fluid.\")\n", - "system.integrator.run(1000)\n", - "logging.info(\"LB fluid warming finished.\")\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", - "\n", - "logging.info(\"Sampling started.\")\n", - "for i in range(LOOPS):\n", - " system.integrator.run(STEPS)\n", - "\n", - "logging.info(\"Sampling finished.\")\n", - "\n", - "correlator.finalize()\n", - "corrdata = correlator.result()\n", - "\n", - "tau = correlator.lag_times()\n", - "msd = corrdata[:, 0] + corrdata[:, 1] + corrdata[:, 2]\n", - "\n", - "rh = system.analysis.calc_rh(chain_start=0, number_of_chains=1, chain_length=N_MONOMERS)[0]" + "N_MONOMERS = np.array([6, 8, 10])\n", + "\n", + "tau_results = []\n", + "msd_results = []\n", + "rh_results = []\n", + "re_results = []\n", + "rg_results = []\n", + "for index, N in enumerate(N_MONOMERS):\n", + " logging.info(\"Polymer size: {}\".format(N))\n", + " # create a linear polymer with Fene bonds\n", + " positions = espressomd.polymer.linear_polymer_positions(n_polymers=1,\n", + " beads_per_chain=N,\n", + " bond_length=1, seed=42,\n", + " min_distance=0.9)\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", + "\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", + " lbf = espressomd.lb.LBFluidGPU(\n", + " kT=1,\n", + " seed=42,\n", + " agrid=1,\n", + " dens=1,\n", + " visc=5,\n", + " tau=TIME_STEP)\n", + " system.actors.add(lbf)\n", + " system.thermostat.set_lb(LB_fluid=lbf, seed=42, gamma=5)\n", + "\n", + " logging.info(\"Warming up the system with LB fluid.\")\n", + " system.integrator.run(1000)\n", + " logging.info(\"Warming up the system with LB fluid finished.\")\n", + "\n", + " # configure correlator\n", + " com_pos = espressomd.observables.ComPosition(ids=range(N))\n", + " correlator = espressomd.accumulators.Correlator(\n", + " 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", + "\n", + " logging.info(\"Sampling started.\")\n", + " rhs = np.zeros(LOOPS)\n", + " res = 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", + " res[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", + "\n", + " logging.info(\"Sampling finished.\")\n", + "\n", + " # store results\n", + " correlator.finalize()\n", + " corrdata = correlator.result()\n", + " tau = correlator.lag_times()\n", + " msd = np.sum(corrdata, axis=1)\n", + " tau_results.append(tau)\n", + " msd_results.append(msd)\n", + " rh_results.append([np.average(rhs), np.std(rhs)])\n", + " re_results.append([np.average(res), np.std(res)])\n", + " rg_results.append([np.average(rgs), np.std(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", + "re_results = np.array(re_results)\n", + "rg_results = np.array(rg_results)\n", + "tau_results = np.array(tau_results)\n", + "msd_results = np.reshape(msd_results, [len(N_MONOMERS),-1])" ] }, { 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 1bc30ae8df9..00000000000 --- a/doc/tutorials/04-lattice_boltzmann/scripts/04-lattice_boltzmann_part3_solution.py +++ /dev/null @@ -1,126 +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 -import espressomd.minimize_energy - -# 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 - espressomd.minimize_energy.steepest_descent( - system, - f_max=1.0, - gamma=10, - max_steps=2000, - max_displacement=0.01) - 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/testsuite/scripts/tutorials/CMakeLists.txt b/testsuite/scripts/tutorials/CMakeLists.txt index 0469cbab9ff..23c6f2432e7 100644 --- a/testsuite/scripts/tutorials/CMakeLists.txt +++ b/testsuite/scripts/tutorials/CMakeLists.txt @@ -35,7 +35,6 @@ 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_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_part3.py b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py index f7ad4624584..36a7608eecb 100644 --- a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py +++ b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py @@ -20,7 +20,7 @@ tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( "@TUTORIALS_DIR@/04-lattice_boltzmann/04-lattice_boltzmann_part3.py", - gpu=True) + LOOPS=500, gpu=True) @skipIfMissingFeatures 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() From 20882e4c7159ce46bcf769cb77813469a491a2cc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Sun, 20 Sep 2020 15:31:28 +0200 Subject: [PATCH 08/16] tutorial 04: part 3: Provide solutions --- .../04-lattice_boltzmann_part3.ipynb | 225 ++++++++++++++++++ .../test_04-lattice_boltzmann_part3.py | 13 +- 2 files changed, 237 insertions(+), 1 deletion(-) 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 deadc78d93b..e5b76b41b40 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb @@ -199,6 +199,231 @@ "msd_results = np.reshape(msd_results, [len(N_MONOMERS),-1])" ] }, + { + "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": "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": [ + "re_exponent, re_prefactor = np.polyfit(np.log(N_MONOMERS), np.log(re_results[:,0]), 1)\n", + "re_prefactor = np.exp(re_prefactor)\n", + "\n", + "fig = plt.figure(figsize=(10, 8))\n", + "x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)\n", + "plt.plot(x, re_prefactor * x**re_exponent, '-',\n", + " label=f'$R_F^{{\\\\mathrm{{fit}}}} = {re_prefactor:.2f} N^{{{re_exponent:.2f}}}$')\n", + "# add error bars assuming uncorrelated samples\n", + "plt.errorbar(N_MONOMERS, re_results[:,0],\n", + " yerr=1.96 * re_results[:,1] / np.sqrt(LOOPS),\n", + " ls='', marker='o', capsize=5, capthick=1,\n", + " label='$R_F^{\\\\mathrm{simulation}}$')\n", + "plt.xlabel('Number of monomers $N$')\n", + "plt.ylabel('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_exponent, rg_prefactor = np.polyfit(np.log(N_MONOMERS), np.log(rg_results[:,0]), 1)\n", + "rg_prefactor = np.exp(rg_prefactor)\n", + "\n", + "fig = plt.figure(figsize=(10, 8))\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=f'$R_g^{{\\\\mathrm{{fit}}}} = {rg_prefactor:.2f} N^{{{rg_exponent:.2f}}}$')\n", + "# add error bars assuming uncorrelated samples\n", + "plt.errorbar(N_MONOMERS, rg_results[:,0],\n", + " yerr=1.96 * rg_results[:,1] / np.sqrt(LOOPS),\n", + " ls='', marker='o', capsize=5, capthick=1,\n", + " label='$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": "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_exponent, rh_prefactor = np.polyfit(np.log(N_MONOMERS), np.log(rh_results[:,0]), 1)\n", + "rh_prefactor = np.exp(rh_prefactor)\n", + "\n", + "fig = plt.figure(figsize=(10, 8))\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=f'$R_h^{{\\\\mathrm{{fit}}}} = {rh_prefactor:.2f} N^{{{rh_exponent:.2f}}}$')\n", + "# add error bars assuming uncorrelated samples\n", + "plt.errorbar(N_MONOMERS, rh_results[:,0],\n", + " yerr=1.96 * rh_results[:,1] / np.sqrt(LOOPS),\n", + " ls='', marker='o', capsize=5, capthick=1,\n", + " label='$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": [ + "Calculate the diffusion coefficient of the polymers." + ] + }, + { + "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(tau_results, msd_results)):\n", + " plt.loglog(tau[1:120], msd[1:120], label=r\"$N=${}\".format(N_MONOMERS[index]))\n", + "plt.loglog(2 * [tau[tau_f_index]], [0, np.max(msd_results)], '-', color='black')\n", + "plt.text(tau[tau_f_index], np.max(msd_results), r'$\\tau_{f}$')\n", + "plt.loglog(2 * [tau[tau_max_index]], [0, np.max(msd_results)], '-', color='black')\n", + "plt.text(tau[tau_max_index], np.max(msd_results), r'$\\tau_{max}$')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "diffusion_results = 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(tau_results, 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_results[index] = a / 6\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Recalling the formula for the diffusion coefficient of a short polymer in the Kirkwood–Zimm model:\n", + "\n", + "$$D = \\frac{D_0}{N} + \\frac{k_B T}{6 \\pi \\eta} \\left\\langle \\frac{1}{R_h} \\right\\rangle$$\n", + "\n", + "where $D_0$ is the monomer diffusion coefficient and $\\eta$ is the viscosity of the fluid." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import scipy.optimize\n", + "\n", + "def kirkwood_zimm(x, a, b, exponent):\n", + " return a / x + b / x**exponent\n", + "\n", + "(a, b), _ = scipy.optimize.curve_fit(\n", + " lambda x, a, b: kirkwood_zimm(x, a, b, rh_exponent),\n", + " N_MONOMERS, diffusion_results)\n", + "\n", + "label = f'''\\\n", + "$D^{{\\\\mathrm{{fit}}}} = \\\n", + " \\\\frac{{{a:.2f}}}{{N}} + \\\n", + " \\\\frac{{{b * 6 * np.pi:.2f} }}{{6\\\\pi}} \\\\cdot \\\n", + " \\\\frac{{{1}}}{{N^{{{rh_exponent:.2f}}}}}$ \\\n", + "'''\n", + "\n", + "fig = plt.figure(figsize=(10, 8))\n", + "x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)\n", + "plt.plot(x, kirkwood_zimm(x, a, b, rh_exponent), '-', label=label)\n", + "plt.plot(N_MONOMERS, diffusion_results, 'o', label='$D^{\\\\mathrm{simulation}}$')\n", + "plt.xlabel('Number of monomers $N$')\n", + "plt.ylabel('Diffusion coefficient')\n", + "plt.legend()\n", + "plt.show()" + ] + }, { "cell_type": "markdown", "metadata": {}, diff --git a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py index 36a7608eecb..62abe1b1a6d 100644 --- a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py +++ b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py @@ -20,13 +20,24 @@ tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( "@TUTORIALS_DIR@/04-lattice_boltzmann/04-lattice_boltzmann_part3.py", - LOOPS=500, gpu=True) + LOOPS=2000, gpu=True) @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.re_exponent, 0.50, msg=msg) + self.assertLess(tutorial.re_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) + if __name__ == "__main__": ut.main() From 15956f91a1509c579df002efa9ef81190dd57578 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Sun, 4 Oct 2020 02:01:41 +0200 Subject: [PATCH 09/16] tutorial 04: part 3: Proper error analysis Integrate the autocorrelation function to remove the variance bias in time series. --- .../04-lattice_boltzmann_part3.ipynb | 102 ++++++++++++++---- .../test_04-lattice_boltzmann_part3.py | 4 +- 2 files changed, 82 insertions(+), 24 deletions(-) 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 e5b76b41b40..76d05106b4f 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb @@ -100,7 +100,7 @@ "tau_results = []\n", "msd_results = []\n", "rh_results = []\n", - "re_results = []\n", + "rf_results = []\n", "rg_results = []\n", "for index, N in enumerate(N_MONOMERS):\n", " logging.info(\"Polymer size: {}\".format(N))\n", @@ -156,7 +156,7 @@ "\n", " logging.info(\"Sampling started.\")\n", " rhs = np.zeros(LOOPS)\n", - " res = 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", @@ -164,7 +164,7 @@ " chain_start=0,\n", " number_of_chains=1,\n", " chain_length=N)[0]\n", - " res[i] = system.analysis.calc_re(\n", + " rfs[i] = system.analysis.calc_re(\n", " chain_start=0,\n", " number_of_chains=1,\n", " chain_length=N)[0]\n", @@ -182,9 +182,9 @@ " msd = np.sum(corrdata, axis=1)\n", " tau_results.append(tau)\n", " msd_results.append(msd)\n", - " rh_results.append([np.average(rhs), np.std(rhs)])\n", - " re_results.append([np.average(res), np.std(res)])\n", - " rg_results.append([np.average(rgs), np.std(rgs)])\n", + " rh_results.append(rhs)\n", + " rf_results.append(rfs)\n", + " rg_results.append(rgs)\n", "\n", " # reset system\n", " system.part.clear()\n", @@ -193,7 +193,7 @@ " system.auto_update_accumulators.clear()\n", "\n", "rh_results = np.array(rh_results)\n", - "re_results = np.array(re_results)\n", + "rf_results = np.array(rf_results)\n", "rg_results = np.array(rg_results)\n", "tau_results = np.array(tau_results)\n", "msd_results = np.reshape(msd_results, [len(N_MONOMERS),-1])" @@ -211,6 +211,64 @@ "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. Due to the short simulation length, it is not possible\n", + " to fit an exponential to the long-time tail. Instead, return a\n", + " 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=f'$\\\\int {variable_label}$ for N={N}')\n", + " summary.append((np.mean(signal), asymptote))\n", + " plt.xlabel('Lag time $\\\\tau / \\\\Delta t$')\n", + " plt.ylabel(f'$\\\\int_{{-\\\\tau}}^{{+\\\\tau}} {variable_label}$')\n", + " plt.legend()\n", + " plt.show()\n", + " return np.array(summary)" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -231,16 +289,16 @@ "metadata": {}, "outputs": [], "source": [ - "re_exponent, re_prefactor = np.polyfit(np.log(N_MONOMERS), np.log(re_results[:,0]), 1)\n", - "re_prefactor = np.exp(re_prefactor)\n", + "rf_summary = standard_error_mean_autocorrelation(rf_results, '\\\\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, 8))\n", "x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)\n", - "plt.plot(x, re_prefactor * x**re_exponent, '-',\n", - " label=f'$R_F^{{\\\\mathrm{{fit}}}} = {re_prefactor:.2f} N^{{{re_exponent:.2f}}}$')\n", - "# add error bars assuming uncorrelated samples\n", - "plt.errorbar(N_MONOMERS, re_results[:,0],\n", - " yerr=1.96 * re_results[:,1] / np.sqrt(LOOPS),\n", + "plt.plot(x, rf_prefactor * x**rf_exponent, '-',\n", + " label=f'$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_F^{\\\\mathrm{simulation}}$')\n", "plt.xlabel('Number of monomers $N$')\n", @@ -269,16 +327,16 @@ "metadata": {}, "outputs": [], "source": [ - "rg_exponent, rg_prefactor = np.polyfit(np.log(N_MONOMERS), np.log(rg_results[:,0]), 1)\n", + "rg_summary = standard_error_mean_autocorrelation(rg_results, '\\\\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, 8))\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=f'$R_g^{{\\\\mathrm{{fit}}}} = {rg_prefactor:.2f} N^{{{rg_exponent:.2f}}}$')\n", - "# add error bars assuming uncorrelated samples\n", - "plt.errorbar(N_MONOMERS, rg_results[:,0],\n", - " yerr=1.96 * rg_results[:,1] / np.sqrt(LOOPS),\n", + "plt.errorbar(N_MONOMERS, rg_summary[:,0],\n", + " yerr=rg_summary[:,1],\n", " ls='', marker='o', capsize=5, capthick=1,\n", " label='$R_g^{\\\\mathrm{simulation}}$')\n", "plt.xlabel('Number of monomers $N$')\n", @@ -311,16 +369,16 @@ }, "outputs": [], "source": [ - "rh_exponent, rh_prefactor = np.polyfit(np.log(N_MONOMERS), np.log(rh_results[:,0]), 1)\n", + "rh_summary = standard_error_mean_autocorrelation(rh_results, '\\\\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, 8))\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=f'$R_h^{{\\\\mathrm{{fit}}}} = {rh_prefactor:.2f} N^{{{rh_exponent:.2f}}}$')\n", - "# add error bars assuming uncorrelated samples\n", - "plt.errorbar(N_MONOMERS, rh_results[:,0],\n", - " yerr=1.96 * rh_results[:,1] / np.sqrt(LOOPS),\n", + "plt.errorbar(N_MONOMERS, rh_summary[:,0],\n", + " yerr=rh_summary[:,1],\n", " ls='', marker='o', capsize=5, capthick=1,\n", " label='$R_h^{\\\\mathrm{simulation}}$')\n", "plt.xlabel('Number of monomers $N$')\n", diff --git a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py index 62abe1b1a6d..5db755402c9 100644 --- a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py +++ b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py @@ -29,8 +29,8 @@ class Tutorial(ut.TestCase): def test_exponents(self): msg = 'The R_F exponent should be close to 0.588' - self.assertGreater(tutorial.re_exponent, 0.50, msg=msg) - self.assertLess(tutorial.re_exponent, 0.85, msg=msg) + 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) From 7d43fd6bf6cb2a420fb43c4586a7d5c655133db5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 6 Oct 2020 17:01:59 +0200 Subject: [PATCH 10/16] tutorial 04: part 3: Green-Kubo diffusion coefficient --- .../04-lattice_boltzmann_part2.ipynb | 2 +- .../04-lattice_boltzmann_part3.ipynb | 247 +++++++++++++++--- .../test_04-lattice_boltzmann_part3.py | 6 + 3 files changed, 215 insertions(+), 40 deletions(-) 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 054918498e6..bf0c970872f 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb @@ -307,7 +307,7 @@ "source": [ "plt.figure(figsize=(10, 8))\n", "plt.xlabel(r'$\\gamma$')\n", - "plt.ylabel('Diffusion coefficient')\n", + "plt.ylabel('Diffusion coefficient [$\\sigma^2/t$]')\n", "gammas = np.linspace(0.9 * min(lb_gammas), 1.1 * max(lb_gammas), 50)\n", "plt.plot(gammas, lbf.kT / gammas, '-', label=r'$k_BT\\gamma^{-1}$')\n", "plt.plot(lb_gammas, diffusion_results, 'o', label='D')\n", 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 76d05106b4f..9d80eba15e9 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb @@ -48,14 +48,15 @@ "\n", "which would create positions for a single polymer with 10 monomers. Please check the documentation for a more detailed description.\n", "\n", - "Furthermore we want to compute the diffusion constant of the polymer for different numbers of monomers.\n", - "For this purpose we can again use the multiple tau correlator. The following script computes the mean\n", - "squared displacement for the center of mass of the polymer as well as the average hydrodynamic radius\n", - "$R_h$, end-to-end distance $R_F$ and radius of gyration $R_g$.\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.\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$." + "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", + "For this purpose we can again use the multiple tau correlator." ] }, { @@ -97,8 +98,10 @@ "\n", "N_MONOMERS = np.array([6, 8, 10])\n", "\n", - "tau_results = []\n", - "msd_results = []\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", @@ -147,12 +150,19 @@ " system.integrator.run(1000)\n", " logging.info(\"Warming up the system with LB fluid finished.\")\n", "\n", - " # configure correlator\n", + " # configure MSD correlator\n", " com_pos = espressomd.observables.ComPosition(ids=range(N))\n", - " correlator = espressomd.accumulators.Correlator(\n", + " com_pos_cor = espressomd.accumulators.Correlator(\n", " 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", + " system.auto_update_accumulators.add(com_pos_cor)\n", + "\n", + " # configure Green-Kubo correlator\n", + " com_vel = espressomd.observables.ComVelocity(ids=range(N))\n", + " com_vel_cor = espressomd.accumulators.Correlator(\n", + " obs1=com_vel, tau_lin=16, tau_max=LOOPS * STEPS, delta_N=10,\n", + " corr_operation=\"scalar_product\", compress1=\"discard1\")\n", + " system.auto_update_accumulators.add(com_vel_cor)\n", "\n", " logging.info(\"Sampling started.\")\n", " rhs = np.zeros(LOOPS)\n", @@ -172,16 +182,15 @@ " chain_start=0,\n", " number_of_chains=1,\n", " chain_length=N)[0]\n", - "\n", " logging.info(\"Sampling finished.\")\n", "\n", " # store results\n", - " correlator.finalize()\n", - " corrdata = correlator.result()\n", - " tau = correlator.lag_times()\n", - " msd = np.sum(corrdata, axis=1)\n", - " tau_results.append(tau)\n", - " msd_results.append(msd)\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", @@ -195,8 +204,18 @@ "rh_results = np.array(rh_results)\n", "rf_results = np.array(rf_results)\n", "rg_results = np.array(rg_results)\n", - "tau_results = np.array(tau_results)\n", - "msd_results = np.reshape(msd_results, [len(N_MONOMERS),-1])" + "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": [ + "We will calculate the means of time series with error bars obtained from\n", + "the correlation-corrected standard error of the mean [8,9]." ] }, { @@ -237,9 +256,10 @@ " '''\n", " Calculate the mean and the correlation-corrected standard error\n", " of the mean of time series by integrating the autocorrelation\n", - " function. Due to the short simulation length, it is not possible\n", - " to fit an exponential to the long-time tail. Instead, return a\n", - " percentile.\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", @@ -269,6 +289,17 @@ " return np.array(summary)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 5.2.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": {}, @@ -391,7 +422,15 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Calculate the diffusion coefficient of the polymers." + "#### 5.2.2 Diffusion coefficient using the MSD method\n", + "\n", + "Calculate the diffusion coefficient of the polymers using the mean-squared displacement.\n", + "\n", + "Recalling that 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." ] }, { @@ -408,12 +447,12 @@ "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", + "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=r\"$N=${}\".format(N_MONOMERS[index]))\n", - "plt.loglog(2 * [tau[tau_f_index]], [0, np.max(msd_results)], '-', color='black')\n", - "plt.text(tau[tau_f_index], np.max(msd_results), r'$\\tau_{f}$')\n", - "plt.loglog(2 * [tau[tau_max_index]], [0, np.max(msd_results)], '-', color='black')\n", - "plt.text(tau[tau_max_index], np.max(msd_results), r'$\\tau_{max}$')\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()" ] @@ -424,17 +463,17 @@ "metadata": {}, "outputs": [], "source": [ - "diffusion_results = np.zeros(len(N_MONOMERS))\n", + "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(tau_results, msd_results)):\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_results[index] = a / 6\n", + " diffusion_msd[index] = a / 6\n", "plt.legend()\n", "plt.show()" ] @@ -443,11 +482,16 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "Plot the dependence of the diffusion coefficient on the hydrodynamic radius.\n", + "\n", "Recalling the formula for the diffusion coefficient of a short polymer in the Kirkwood–Zimm model:\n", "\n", "$$D = \\frac{D_0}{N} + \\frac{k_B T}{6 \\pi \\eta} \\left\\langle \\frac{1}{R_h} \\right\\rangle$$\n", "\n", - "where $D_0$ is the monomer diffusion coefficient and $\\eta$ is the viscosity of the fluid." + "where $D_0$ is the monomer diffusion coefficient and $\\eta$ is the viscosity of the fluid.\n", + "\n", + "Hint: use $D = \\alpha_1 N^{-1} + \\alpha_2 N^{-\\beta}$ with `rh_exponent` for $\\beta$\n", + "and solve for $\\alpha_1, \\alpha_2$." ] }, { @@ -463,25 +507,148 @@ "\n", "(a, b), _ = scipy.optimize.curve_fit(\n", " lambda x, a, b: kirkwood_zimm(x, a, b, rh_exponent),\n", - " N_MONOMERS, diffusion_results)\n", + " N_MONOMERS, diffusion_msd)\n", "\n", "label = f'''\\\n", "$D^{{\\\\mathrm{{fit}}}} = \\\n", " \\\\frac{{{a:.2f}}}{{N}} + \\\n", - " \\\\frac{{{b * 6 * np.pi:.2f} }}{{6\\\\pi}} \\\\cdot \\\n", + " \\\\frac{{{b * 6 * np.pi:.3f} }}{{6\\\\pi}} \\\\cdot \\\n", + " \\\\frac{{{1}}}{{N^{{{rh_exponent:.2f}}}}}$ \\\n", + "'''\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, kirkwood_zimm(x, a, b, rh_exponent), '-', label=label)\n", + "plt.plot(N_MONOMERS, diffusion_msd, 'o', label='$D^{\\\\mathrm{simulation}}$')\n", + "plt.xlabel('Number of monomers $N$')\n", + "plt.ylabel('Diffusion coefficient [$\\\\sigma^2/t$]')\n", + "plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 4))\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 5.2.3 Diffusion coefficient using the Green–Kubo method\n", + "\n", + "Plot the autocorrelation function and check that the decay is roughly exponential.\n", + "\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=f'$R(\\\\tau)$ for $N = {N}$')\n", + "plt.xlabel('$\\\\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", + "$$D = \\frac{1}{3} \\int_0^{+\\infty} \\left<\\vec{v_c}(\\tau)\\cdot\\vec{v_c}(0)\\right>\\, \\mathrm{d}\\tau$$\n", + "\n", + "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=f'$D(\\\\tau_{{\\\\mathrm{{int}}}})$ for $N = {N}$')\n", + " diffusion_gk.append(np.mean(y[10:]))\n", + "plt.xlabel('$\\\\tau_{\\\\mathrm{int}}$')\n", + "plt.ylabel('$\\\\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=(1, 4))\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": [ + "(a, b), _ = scipy.optimize.curve_fit(\n", + " lambda x, a, b: kirkwood_zimm(x, a, b, rh_exponent),\n", + " N_MONOMERS, diffusion_gk)\n", + "\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", "\n", "fig = plt.figure(figsize=(10, 8))\n", "x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)\n", "plt.plot(x, kirkwood_zimm(x, a, b, rh_exponent), '-', label=label)\n", - "plt.plot(N_MONOMERS, diffusion_results, 'o', label='$D^{\\\\mathrm{simulation}}$')\n", + "plt.plot(N_MONOMERS, diffusion_gk, 'o', label='$D^{\\\\mathrm{simulation}}$')\n", "plt.xlabel('Number of monomers $N$')\n", - "plt.ylabel('Diffusion coefficient')\n", + "plt.ylabel('Diffusion coefficient [$\\\\sigma^2/t$]')\n", + "plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 4))\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}%')" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -494,7 +661,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/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py index 5db755402c9..1d57a5dc095 100644 --- a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py +++ b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py @@ -17,6 +17,7 @@ import unittest as ut import importlib_wrapper +import numpy as np tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( "@TUTORIALS_DIR@/04-lattice_boltzmann/04-lattice_boltzmann_part3.py", @@ -38,6 +39,11 @@ def test_exponents(self): 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.1) + np.testing.assert_allclose(tutorial.diffusion_gk, reference, rtol=0.1) + if __name__ == "__main__": ut.main() From 02e102fd3ae55a25f3eca6db030e136b11294abf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 6 Oct 2020 18:01:17 +0200 Subject: [PATCH 11/16] tutorial 04: Add tutor notes --- .../04-lattice_boltzmann/tutor_notes.md | 25 +++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 doc/tutorials/04-lattice_boltzmann/tutor_notes.md 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..1dbbca2073c --- /dev/null +++ b/doc/tutorials/04-lattice_boltzmann/tutor_notes.md @@ -0,0 +1,25 @@ +# Part 2: Brownian motion of a sphere in a fluid + +## Physical learning goals + +* Learning the basics of Brownian motion +* Getting a short introduction to the Flory theory of polymers +* How to calculate the diffusion coefficient of a sphere + +## Espresso learning goals + +* Set up a LB fluid with particle coupling + +# Part 3: Brownian motion of a polymer in a fluid + +## Physical learning goals + +* Model the Brownian motion of a polymer +* Learn about the hydrodynamic radius, radius of gyration, end-to-end distance +* Reproduce 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 + +## Espresso learning goals + +* Using correlators From f7cd6abede7772286873df0d8fdd402c0b9f2ffe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 7 Oct 2020 14:22:00 +0200 Subject: [PATCH 12/16] tutorial 04: part 4: Switch exercise The participants should focus on setting up an LB fluid instead of writing the matplotlib code. --- .../04-lattice_boltzmann_part4.ipynb | 172 +++++++++++++----- .../04-lattice_boltzmann/tutor_notes.md | 10 + .../test_04-lattice_boltzmann_part4.py | 2 +- 3 files changed, 141 insertions(+), 43 deletions(-) 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..4bf2132d5d8 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part4.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part4.ipynb @@ -41,6 +41,13 @@ "applied to every node." ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### 6.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": [ + "#### 6.2 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,87 @@ }, "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": [ + "#### 6.3 Simulation and data analysis\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": [ + "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/tutor_notes.md b/doc/tutorials/04-lattice_boltzmann/tutor_notes.md index 1dbbca2073c..7e6c3ff47eb 100644 --- a/doc/tutorials/04-lattice_boltzmann/tutor_notes.md +++ b/doc/tutorials/04-lattice_boltzmann/tutor_notes.md @@ -23,3 +23,13 @@ ## Espresso learning goals * Using correlators + +# 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/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) From 59984d130fcb23d336e33f48e31f7fa1fcab11ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 7 Oct 2020 15:36:04 +0200 Subject: [PATCH 13/16] tutorial 04: Create solution cells Also give more details on tasks, add links to the user guide and renumber sections. --- .../04-lattice_boltzmann_part2.ipynb | 64 ++++-- .../04-lattice_boltzmann_part3.ipynb | 203 +++++++++++++----- .../04-lattice_boltzmann_part4.ipynb | 10 +- .../04-lattice_boltzmann/tutor_notes.md | 17 +- 4 files changed, 216 insertions(+), 78 deletions(-) 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 bf0c970872f..d8e1b81f484 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb @@ -11,7 +11,7 @@ "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", @@ -71,14 +71,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": [ + "### 2. Simulating the Brownian motion" ] }, { "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:" + "We will simulate the diffusion of a single particle that is coupled to an LB fluid by the point coupling method:" ] }, { @@ -111,9 +152,6 @@ "\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 coefficients\n", "lb_gammas = [1.0, 2.0, 4.0, 10.0]\n", "tau_results = []\n", @@ -123,14 +161,13 @@ " 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", + "\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(\n", - " obs1=pos, tau_lin=16, tau_max=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", @@ -146,7 +183,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### Plotting the results" + "### 3. Data analysis\n", + "#### 3.1 Plotting the results" ] }, { @@ -167,7 +205,7 @@ " # 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", " # the high variance for larger lag times.\n", - " plt.loglog(tau[1:100], msd[1:100], label=r\"$\\gamma=${:.1f}\".format(lb_gammas[index]))\n", + " plt.loglog(tau[1:100], msd[1:100], label=r'$\\gamma=${:.1f}'.format(lb_gammas[index]))\n", "plt.legend()\n", "plt.show()" ] @@ -176,7 +214,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### Calculating the diffusion coefficient\n", + "#### 3.2 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", "The mean squared displacement is calculated during the simulation via a multiple-tau\n", 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 9d80eba15e9..41b0ff5ae2e 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb @@ -11,7 +11,7 @@ "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", @@ -46,17 +46,114 @@ " 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.\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", "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", - "For this purpose we can again use the multiple tau correlator." + "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": [ + "### 2. Simulating the polymer" ] }, { @@ -83,6 +180,7 @@ "TIME_STEP = 0.01\n", "LOOPS = 4000\n", "STEPS = 100\n", + "POLYMER_PARAMS = {'n_polymers': 1, 'bond_length': 1, 'seed': 42, 'min_distance': 0.9}\n", "\n", "# System setup\n", "system = espressomd.System(box_l=[12.0, 12.0, 12.0])\n", @@ -107,16 +205,7 @@ "rg_results = []\n", "for index, N in enumerate(N_MONOMERS):\n", " logging.info(\"Polymer size: {}\".format(N))\n", - " # create a linear polymer with Fene bonds\n", - " positions = espressomd.polymer.linear_polymer_positions(n_polymers=1,\n", - " beads_per_chain=N,\n", - " bond_length=1, seed=42,\n", - " min_distance=0.9)\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", + " build_polymer(system, N, POLYMER_PARAMS, fene)\n", "\n", " logging.info(\"Warming up the polymer chain.\")\n", " system.time_step = 0.002\n", @@ -151,17 +240,11 @@ " logging.info(\"Warming up the system with LB fluid finished.\")\n", "\n", " # configure MSD correlator\n", - " com_pos = espressomd.observables.ComPosition(ids=range(N))\n", - " com_pos_cor = espressomd.accumulators.Correlator(\n", - " obs1=com_pos, tau_lin=16, tau_max=LOOPS * STEPS, delta_N=1,\n", - " corr_operation=\"square_distance_componentwise\", compress1=\"discard1\")\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 = espressomd.observables.ComVelocity(ids=range(N))\n", - " com_vel_cor = espressomd.accumulators.Correlator(\n", - " obs1=com_vel, tau_lin=16, tau_max=LOOPS * STEPS, delta_N=10,\n", - " corr_operation=\"scalar_product\", compress1=\"discard1\")\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", @@ -214,6 +297,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "### 3. Data analysis\n", + "\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]." ] @@ -280,10 +365,11 @@ " 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=f'$\\\\int {variable_label}$ for N={N}')\n", + " label=rf'$\\int {variable_label}$ for N={N}')\n", " summary.append((np.mean(signal), asymptote))\n", - " plt.xlabel('Lag time $\\\\tau / \\\\Delta t$')\n", - " plt.ylabel(f'$\\\\int_{{-\\\\tau}}^{{+\\\\tau}} {variable_label}$')\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)" @@ -293,7 +379,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### 5.2.1 Distance-based macromolecular properties\n", + "#### 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", @@ -320,20 +406,20 @@ "metadata": {}, "outputs": [], "source": [ - "rf_summary = standard_error_mean_autocorrelation(rf_results, '\\\\operatorname{acf}(R_F)')\n", + "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, 8))\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=f'$R_F^{{\\\\mathrm{{fit}}}} = {rf_prefactor:.2f} N^{{{rf_exponent:.2f}}}$')\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_F^{\\\\mathrm{simulation}}$')\n", + " label=r'$R_F^{\\mathrm{simulation}}$')\n", "plt.xlabel('Number of monomers $N$')\n", - "plt.ylabel('End-to-end distance [$\\sigma$]')\n", + "plt.ylabel(r'End-to-end distance [$\\sigma$]')\n", "plt.legend()\n", "plt.show()" ] @@ -358,24 +444,33 @@ "metadata": {}, "outputs": [], "source": [ - "rg_summary = standard_error_mean_autocorrelation(rg_results, '\\\\operatorname{acf}(R_g)')\n", + "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, 8))\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=f'$R_g^{{\\\\mathrm{{fit}}}} = {rg_prefactor:.2f} N^{{{rg_exponent:.2f}}}$')\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_g^{\\\\mathrm{simulation}}$')\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": {}, @@ -400,18 +495,18 @@ }, "outputs": [], "source": [ - "rh_summary = standard_error_mean_autocorrelation(rh_results, '\\\\operatorname{acf}(R_h)')\n", + "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, 8))\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=f'$R_h^{{\\\\mathrm{{fit}}}} = {rh_prefactor:.2f} N^{{{rh_exponent:.2f}}}$')\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_h^{\\\\mathrm{simulation}}$')\n", + " label=r'$R_h^{\\mathrm{simulation}}$')\n", "plt.xlabel('Number of monomers $N$')\n", "plt.ylabel('Hydrodynamic radius [$\\sigma$]')\n", "plt.legend()\n", @@ -422,7 +517,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### 5.2.2 Diffusion coefficient using the MSD method\n", + "#### 3.2 Diffusion coefficient using the MSD method\n", "\n", "Calculate the diffusion coefficient of the polymers using the mean-squared displacement.\n", "\n", @@ -448,7 +543,7 @@ "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=r\"$N=${}\".format(N_MONOMERS[index]))\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", @@ -519,10 +614,10 @@ "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, kirkwood_zimm(x, a, b, rh_exponent), '-', label=label)\n", - "plt.plot(N_MONOMERS, diffusion_msd, 'o', label='$D^{\\\\mathrm{simulation}}$')\n", + "plt.plot(N_MONOMERS, diffusion_msd, 'o', label=r'$D^{\\mathrm{simulation}}$')\n", "plt.xlabel('Number of monomers $N$')\n", - "plt.ylabel('Diffusion coefficient [$\\\\sigma^2/t$]')\n", - "plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 4))\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()" ] @@ -531,7 +626,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### 5.2.3 Diffusion coefficient using the Green–Kubo method\n", + "#### 3.3 Diffusion coefficient using the Green–Kubo method\n", "\n", "Plot the autocorrelation function and check that the decay is roughly exponential.\n", "\n", @@ -555,8 +650,8 @@ " 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=f'$R(\\\\tau)$ for $N = {N}$')\n", - "plt.xlabel('$\\\\tau$')\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()" @@ -587,11 +682,11 @@ "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=f'$D(\\\\tau_{{\\\\mathrm{{int}}}})$ for $N = {N}$')\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('$\\\\tau_{\\\\mathrm{int}}$')\n", - "plt.ylabel('$\\\\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=(1, 4))\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()" ] @@ -623,10 +718,10 @@ "fig = plt.figure(figsize=(10, 8))\n", "x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)\n", "plt.plot(x, kirkwood_zimm(x, a, b, rh_exponent), '-', label=label)\n", - "plt.plot(N_MONOMERS, diffusion_gk, 'o', label='$D^{\\\\mathrm{simulation}}$')\n", + "plt.plot(N_MONOMERS, diffusion_gk, 'o', label=r'$D^{\\mathrm{simulation}}$')\n", "plt.xlabel('Number of monomers $N$')\n", - "plt.ylabel('Diffusion coefficient [$\\\\sigma^2/t$]')\n", - "plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 4))\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()" ] 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 4bf2132d5d8..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", @@ -45,7 +45,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### 6.1 setting up the system" + "### 1. Setting up the system" ] }, { @@ -82,7 +82,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### 6.2 Setting up the lattice-Boltzmann fluid\n", + "#### 1.1 Setting up the lattice-Boltzmann fluid\n", "\n", "We will now create a lattice-Boltzmann fluid confined between two walls." ] @@ -179,7 +179,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### 6.3 Simulation and data analysis\n", + "### 2. Simulation\n", "\n", "We will now simulate the fluid flow until we reach the steady state." ] @@ -198,6 +198,8 @@ "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." ] }, diff --git a/doc/tutorials/04-lattice_boltzmann/tutor_notes.md b/doc/tutorials/04-lattice_boltzmann/tutor_notes.md index 7e6c3ff47eb..7f63111de9e 100644 --- a/doc/tutorials/04-lattice_boltzmann/tutor_notes.md +++ b/doc/tutorials/04-lattice_boltzmann/tutor_notes.md @@ -2,27 +2,30 @@ ## Physical learning goals -* Learning the basics of Brownian motion * Getting a short introduction to the Flory theory of polymers -* How to calculate the diffusion coefficient of a sphere +* Learning the basics of Brownian motion +* Modeling the Brownian motion of a spherical object +* Calculating the diffusion coefficient of a spherical object ## Espresso learning goals -* Set up a LB fluid with particle coupling +* Setting up a LB fluid with particle coupling # Part 3: Brownian motion of a polymer in a fluid ## Physical learning goals -* Model the Brownian motion of a polymer -* Learn about the hydrodynamic radius, radius of gyration, end-to-end distance -* Reproduce theoretical results of the Flory theory of polymers +* 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 ## Espresso learning goals -* Using correlators +* 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 From b71581d258fe73bee02a6a33d4aa590915ee5117 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 7 Oct 2020 16:24:32 +0200 Subject: [PATCH 14/16] tutorial 04: part 1: Update text, add links --- .../04-lattice_boltzmann_part1.ipynb | 31 ++++++++++--------- 1 file changed, 16 insertions(+), 15 deletions(-) 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 b9c7b6bc2b2..9501412c6ff 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb @@ -29,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", @@ -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**." ] }, { @@ -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", @@ -272,7 +272,7 @@ "+ 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", @@ -283,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", @@ -315,8 +315,9 @@ "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 result 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", From c7d42c89ab9a9907394bdbed62b34e8453e95cb7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 7 Oct 2020 17:45:27 +0200 Subject: [PATCH 15/16] tutorial 04: Polymer diffusion with implicit solvent The choice of solvent implementation only affects the Rouse vs. Zimm regime. Chain properties such as the hydrodynamics radius, end-to-end distance and radius of gyration do not change. --- .../04-lattice_boltzmann_part2.ipynb | 49 ++++---- .../04-lattice_boltzmann_part3.ipynb | 110 +++++++++++------- .../04-lattice_boltzmann/tutor_notes.md | 15 ++- testsuite/scripts/tutorials/CMakeLists.txt | 5 +- .../test_04-lattice_boltzmann_part2.py | 5 +- .../test_04-lattice_boltzmann_part3.py | 12 +- 6 files changed, 122 insertions(+), 74 deletions(-) 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 d8e1b81f484..ffda50b6fe0 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb @@ -13,13 +13,18 @@ "source": [ "## 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", - "\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", + "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\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", @@ -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 polymer in an\n", + "implicit solvent and in an LB fluid, 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." @@ -119,7 +123,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We will simulate the diffusion of a single particle that is coupled to an LB fluid by the point coupling method:" + "We will simulate the diffusion of a single particle that is coupled to an implicit solvent." ] }, { @@ -140,6 +144,7 @@ "logging.basicConfig(level=logging.INFO, stream=sys.stdout)\n", "\n", "# Constants\n", + "KT = 1.1\n", "STEPS = 400000\n", "\n", "# System setup\n", @@ -147,20 +152,17 @@ "system.time_step = 0.01\n", "system.cell_system.skin = 0.4\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", "system.part.add(pos=[0, 0, 0])\n", "\n", "# Run for different friction coefficients\n", - "lb_gammas = [1.0, 2.0, 4.0, 10.0]\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", + " system.thermostat.set_langevin(kT=KT, gamma=gamma, seed=42)\n", "\n", " logging.info(\"Equilibrating the system.\")\n", " system.integrator.run(1000)\n", @@ -205,7 +207,7 @@ " # 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", " # the high variance for larger lag times.\n", - " plt.loglog(tau[1:100], msd[1:100], label=r'$\\gamma=${:.1f}'.format(lb_gammas[index]))\n", + " plt.loglog(tau[1:100], msd[1:100], label=r'$\\gamma=${:.1f}'.format(gammas[index]))\n", "plt.legend()\n", "plt.show()" ] @@ -269,7 +271,7 @@ " 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(lb_gammas[index]))\n", + " label=r'$\\gamma=${:.1f}'.format(gammas[index]))\n", "plt.legend()\n", "plt.show()" ] @@ -316,7 +318,7 @@ " 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(lb_gammas[index]))\n", + " label=r'$\\gamma=${:.1f}'.format(gammas[index]))\n", " diffusion_results.append(a / 6)\n", "\n", "plt.legend()\n", @@ -346,9 +348,10 @@ "plt.figure(figsize=(10, 8))\n", "plt.xlabel(r'$\\gamma$')\n", "plt.ylabel('Diffusion coefficient [$\\sigma^2/t$]')\n", - "gammas = np.linspace(0.9 * min(lb_gammas), 1.1 * max(lb_gammas), 50)\n", - "plt.plot(gammas, lbf.kT / gammas, '-', label=r'$k_BT\\gamma^{-1}$')\n", - "plt.plot(lb_gammas, diffusion_results, 'o', label='D')\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()" ] 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 41b0ff5ae2e..e99c77ca458 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb @@ -149,6 +149,36 @@ "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": {}, @@ -180,7 +210,11 @@ "TIME_STEP = 0.01\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", "\n", "# System setup\n", "system = espressomd.System(box_l=[12.0, 12.0, 12.0])\n", @@ -225,15 +259,10 @@ "\n", " system.thermostat.turn_off()\n", "\n", - " lbf = espressomd.lb.LBFluidGPU(\n", - " kT=1,\n", - " seed=42,\n", - " agrid=1,\n", - " dens=1,\n", - " visc=5,\n", - " tau=TIME_STEP)\n", - " system.actors.add(lbf)\n", - " system.thermostat.set_lb(LB_fluid=lbf, seed=42, gamma=5)\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 LB fluid.\")\n", " system.integrator.run(1000)\n", @@ -584,9 +613,11 @@ "$$D = \\frac{D_0}{N} + \\frac{k_B T}{6 \\pi \\eta} \\left\\langle \\frac{1}{R_h} \\right\\rangle$$\n", "\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", - "Hint: use $D = \\alpha_1 N^{-1} + \\alpha_2 N^{-\\beta}$ with `rh_exponent` for $\\beta$\n", - "and solve for $\\alpha_1, \\alpha_2$." + "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." ] }, { @@ -597,23 +628,35 @@ "source": [ "import scipy.optimize\n", "\n", - "def kirkwood_zimm(x, a, b, exponent):\n", - " return a / x + b / x**exponent\n", - "\n", - "(a, b), _ = scipy.optimize.curve_fit(\n", - " lambda x, a, b: kirkwood_zimm(x, a, b, rh_exponent),\n", - " N_MONOMERS, diffusion_msd)\n", - "\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", + "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 = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)\n", - "plt.plot(x, kirkwood_zimm(x, a, b, rh_exponent), '-', label=label)\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", @@ -704,20 +747,9 @@ "metadata": {}, "outputs": [], "source": [ - "(a, b), _ = scipy.optimize.curve_fit(\n", - " lambda x, a, b: kirkwood_zimm(x, a, b, rh_exponent),\n", - " N_MONOMERS, diffusion_gk)\n", - "\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", - "\n", "fig = plt.figure(figsize=(10, 8))\n", - "x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)\n", - "plt.plot(x, kirkwood_zimm(x, a, b, rh_exponent), '-', label=label)\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", diff --git a/doc/tutorials/04-lattice_boltzmann/tutor_notes.md b/doc/tutorials/04-lattice_boltzmann/tutor_notes.md index 7f63111de9e..ac7a0576080 100644 --- a/doc/tutorials/04-lattice_boltzmann/tutor_notes.md +++ b/doc/tutorials/04-lattice_boltzmann/tutor_notes.md @@ -1,3 +1,13 @@ +# 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 @@ -7,10 +17,6 @@ * Modeling the Brownian motion of a spherical object * Calculating the diffusion coefficient of a spherical object -## Espresso learning goals - -* Setting up a LB fluid with particle coupling - # Part 3: Brownian motion of a polymer in a fluid ## Physical learning goals @@ -20,6 +26,7 @@ * 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 diff --git a/testsuite/scripts/tutorials/CMakeLists.txt b/testsuite/scripts/tutorials/CMakeLists.txt index 23c6f2432e7..ef3259fda00 100644 --- a/testsuite/scripts/tutorials/CMakeLists.txt +++ b/testsuite/scripts/tutorials/CMakeLists.txt @@ -33,8 +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_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 d71cc5fa847..805b5445459 100644 --- a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py +++ b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part2.py @@ -21,8 +21,7 @@ import scipy.optimize tutorial, skipIfMissingFeatures = importlib_wrapper.configure_and_import( - "@TUTORIALS_DIR@/04-lattice_boltzmann/04-lattice_boltzmann_part2.py", - gpu=True, STEPS=8000, cutoff_limit=66) + "@TUTORIALS_DIR@/04-lattice_boltzmann/04-lattice_boltzmann_part2.py") @skipIfMissingFeatures @@ -40,7 +39,7 @@ def test_ballistic_regime(self): def test_diffusion_coefficient(self): D_val = tutorial.diffusion_results - D_ref = tutorial.lbf.kT / np.array(tutorial.lb_gammas) + D_ref = tutorial.KT / np.array(tutorial.gammas) np.testing.assert_allclose(D_val, D_ref, rtol=0, atol=0.1) diff --git a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py index 1d57a5dc095..56075de1c3e 100644 --- a/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py +++ b/testsuite/scripts/tutorials/test_04-lattice_boltzmann_part3.py @@ -19,9 +19,14 @@ 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", - LOOPS=2000, gpu=True) + script_suffix="@TEST_SUFFIX@", **params) @skipIfMissingFeatures @@ -41,8 +46,9 @@ def test_exponents(self): def test_diffusion_coefficients(self): reference = [0.0363, 0.0269, 0.0234] - np.testing.assert_allclose(tutorial.diffusion_msd, reference, rtol=0.1) - np.testing.assert_allclose(tutorial.diffusion_gk, reference, rtol=0.1) + 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__": From 8469df073f28eda549e36362d2455e348fd909c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 7 Oct 2020 20:54:17 +0200 Subject: [PATCH 16/16] tutorial 04: Remove mentions of LB --- .../04-lattice_boltzmann_part2.ipynb | 10 +++++----- .../04-lattice_boltzmann_part3.ipynb | 7 +++++-- 2 files changed, 10 insertions(+), 7 deletions(-) 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 ffda50b6fe0..d23a2e6484f 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part2.ipynb @@ -64,8 +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 found in virtually any simulation textbook, like [7]. We will set up a polymer in an\n", - "implicit solvent and 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." @@ -138,7 +138,6 @@ "\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", @@ -217,8 +216,9 @@ "metadata": {}, "source": [ "#### 3.2 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", + "\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?" 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 e99c77ca458..12a76ef4074 100644 --- a/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb +++ b/doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part3.ipynb @@ -215,6 +215,9 @@ "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=[12.0, 12.0, 12.0])\n", @@ -264,9 +267,9 @@ " elif POLYMER_MODEL == 'Zimm':\n", " solvent_lbm(system, KT, GAMMA)\n", "\n", - " logging.info(\"Warming up the system with LB fluid.\")\n", + " logging.info(\"Warming up the system with the fluid.\")\n", " system.integrator.run(1000)\n", - " logging.info(\"Warming up the system with LB fluid finished.\")\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",