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.

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,27 +144,25 @@
     "logging.basicConfig(level=logging.INFO, stream=sys.stdout)\n",
     "\n",
     "# Constants\n",
+    "KT = 1.1\n",
     "STEPS = 400000\n",
     "\n",
     "# System setup\n",
     "system = espressomd.System(box_l=[16] * 3)\n",
     "system.time_step = 0.01\n",
     "system.cell_system.skin = 0.4\n",
     "\n",
-    "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()"
    ]

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",

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
 

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")

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)
 
 

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__":