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 ν is the Flory exponent
that plays a very prominent role in polymer physics. It has a value of $∼ 3/5$ in good
solvent conditions in 3D. Discussions of these scaling laws can be found in polymer
physics textbooks like [4–6].

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 long enough diffuse more like 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
hydrodynamic radius. It is defined as:

\begin{equation} \langle \frac{1}{R_h} \rangle = \langle \frac{1}{N^2}\sum_{i\neq j} \frac{1}{\left| r_i - r_j \right|} \rangle \end{equation}

This hydrodynamic radius exhibits the scaling law $R_h \propto N^{\nu}$
and the diffusion coefficient of long polymer is proportional to its inverse.

For shorter polymers there is a transition region. It can be described
by the Kirkwood-Zimm model:
\begin{equation}
D=\frac{D_0}{N} + \frac{k_B T}{6 \pi \eta } \langle \frac{1}{R_h} \rangle
\end{equation}

Here $D_0$ is the monomer diffusion coefficient and $\eta$ the viscosity of the fluid. For a finite system size the second part of the diffusion is subject of a $1/L$ finite size effect, because hydrodynamic interactions are proportional to the inverse distance and thus long ranged. It can be taken into account by a correction:

\begin{equation}
D=\frac{D_0}{N} + \frac{k_B T}{6 \pi \eta } \langle \frac{1}{R_h} \rangle \left( 1- \langle\frac{R_h}{L} \rangle \right)
\end{equation}
It is quite difficult to prove this formula with good accuracy. It will
need quite some computer time and a careful analysis. So please don't be
too disappointed if you don't manage to do so.

We want to determine the diffusion coefficient from the mean square
distance that a particle travels in the time $t$. For large $t$ it is
be proportional to the time and the diffusion coefficient occurs as
prefactor:
\begin{equation}
\frac{\partial \langle r^2 \left(t\right)\rangle}{\partial t} = 2 d D.
\end{equation}
Here $d$ denotes the dimensionality of the system, in our case 3.This equation can be
found in virtually any simulation textbook, like [7]. We will therefore set up a polymer
in an LB fluid, simulate for an appropriate amount of time, calculate the mean square
displacement as a function of time and obtain the diffusion coefficient from a linear
fit. However we make a couple of steps in between and divide the full problem into
subproblems that allow to (hopefully) fully understand the process.

Our first step is to investigate the diffusion of a single particle that is coupled to an LB fluid by the point coupling method. Take a look at the script `single_particle_diffusion.py`. The script takes the LB-friction coefficient as an argument. Start with an friction coefficient of 1.0:

In [1]:

```
# The friction coefficient has to be changed inside the code as sys.argv[1]
# does not work in iPython.
from espressomd import System, lb
from espressomd.observables import ParticlePositions
from espressomd.accumulators import Correlator
import numpy as np
import sys
# Constants
loops = 40000
steps = 10
time_step = 0.01
runtime = loops*steps*time_step
box_l = 16
lb_friction = 1.0 # The friction has been set to be 1.0 as default
# Change this value for further simulations
# System setup
system = System(box_l = [box_l] * 3)
system.set_random_state_PRNG()
np.random.seed(seed = system.seed)
system.time_step = time_step
system.cell_system.skin = 0.4
lbf = lb.LBFluidGPU(agrid=1, dens=1, visc=5, tau=0.01, fric=lb_friction)
system.actors.add(lbf)
system.thermostat.set_lb(kT=1)
system.part.add(pos=[0, 0, 0])
## perform a couple of steps to come to equilbrium
print("Equilibrating the system.")
system.integrator.run(1000)
print("Equlibration finished.")
# Setup observable correlator
pos = ParticlePositions(ids=(0,))
c = Correlator(obs1=pos, tau_lin = 16, tau_max = 1000, delta_N = 1,
corr_operation="square_distance_componentwise", compress1="discard1")
system.auto_update_accumulators.add(c)
print("Sampling started.")
for i in range(loops):
system.integrator.run(steps)
if i % 1e2 == 0:
sys.stdout.write("\rSampling: %05i"%i)
sys.stdout.flush()
print("Sampling finished.")
c.finalize()
corrdata = c.result()
corr = np.zeros((corrdata.shape[0],2))
corr[:,0] = corrdata[:,0]
corr[:,1] = (corrdata[:,2] + corrdata[:,3] + corrdata[:,4]) / 3
np.savetxt("./msd_"+str(lb_friction)+".dat", corr)
```

In this script an LB fluid and a single particle are created and thermalized. The
random forces on the particle and within the LB fluid will cause the particle to move.
The mean squared displacement is calculated during the simulation via a multiple-tau
correlator. Run the simulation script and plot the output data `msd_1.0.dat`. To load the file into a numpy array, one can use `numpy.loadtxt`. Zoom in on the origin of the
plot. What do you see for short times? What do you see on a longer time scale? Produce
a double-logarithmic plot to assess the power law.

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 coefficient. Run the simulation again with different values for the friction coefficient, e.g., 1, 2, 4, and 10. Calculate the diffusion coefficient for all cases and plot them as a function of γ. What relation do you observe?

[1] S Succi. *The lattice Boltzmann equation for fluid dynamics and beyond. *Clarendon Press, Oxford, 2001.

[2] B. Dünweg and A. J. C. Ladd. *Advanced Computer Simulation Approaches for Soft Matter Sciences III*, chapter II, pages 89–166. Springer, 2009.

[3] B. Dünweg, U. Schiller, and A.J.C. Ladd. Statistical mechanics of the fluctuating lattice-boltzmann equation. *Phys. Rev. E*, 76:36704, 2007.

[4] P. G. de Gennes. *Scaling Concepts in Polymer Physics*. Cornell University Press, Ithaca, NY, 1979.

[5] M. Doi. *Introduction do Polymer Physics.* Clarendon Press, Oxford, 1996.

[6] Michael Rubinstein and Ralph H. Colby. *Polymer Physics.* Oxford University Press, Oxford, UK, 2003.

[7] Daan Frenkel and Berend Smit. *Understanding Molecular Simulation.* Academic Press, San Diego, second edition, 2002.