**ESPResSo**. For this tutorial, you should have fundamental knowledge of the lattice-Boltzmann method and Langevin dynamics. If you are unfamiliar with those, you can go through the respective tutorials in the `lattice_boltzmann`

and `langevin_dynamics`

folders.

In these exercises we want 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 that plays a very prominent role in polymer physics. It has a value of $\sim 3/5$ in good solvent conditions in 3D. Discussions on these scaling laws can be found in polymer physics textbooks like [1], [2], and [3, chapter 8].

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 hydrodynamic radius. It is defined as:

$$ \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 $$This hydrodynamic radius exhibits the scaling law $R_h \propto N^{\nu}$ and the diffusion coefficient of a long polymer is proportional to its inverse $R_h$.

The diffusion coefficient $D$ of a spherical particle in a carrier fluid experiencing drag is related to the friction coefficient $\zeta$ via the Einstein relation:

$$ D = \frac{k_{\mathrm{B}}T}{\zeta}, $$with $k_{\mathrm{B}}$ the Boltzmann constant and $T$ the temperature. For a sphere of radius $R$ moving in a fluid of viscosity $\eta$, the friction coefficient is obtained via the Stokes law:

$$ \zeta = 6\pi\eta R. $$Combining both equations yields the Stokes–Einstein relation:

$$ D = \frac{k_{\mathrm{B}}T}{6\pi\eta R}. $$The simplest description of a polymer is the Rouse model, where beads are connected by springs. All beads experience a drag from the solvent, and the friction coefficient $\gamma$ is identical for all beads. The solvent flows freely between beads and hydrodynamic interactions are neglected. The diffusion coefficient takes the following form:

$$ D_{\mathrm{R}} = \frac{D_0}{N} = \frac{k_{\mathrm{B}}T}{\gamma N}, $$where $D_0$ is the diffusion coefficient of a single bead.

To account for hydrodynamic interactions mediated by the solvent, i.e. the transport of solvent in contact with the beads and the correlation in the motion of beads due to the carried solvent, the Zimm model was created. For an ideal chain, it takes the following form:

$$ D_{\mathrm{Z}} = \frac{8}{3\sqrt{6\pi^3}}\frac{k_B T}{\eta R} \simeq 0.196\frac{k_B T}{\eta b N^{\nu}}, $$with $R$ the radius of the polymer and $b$ the length of the spring connecting the beads.

For shorter polymers there is a transition region. It can be described by the Kirkwood–Zimm model:

$$ D=\frac{D_0}{N} + \frac{k_B T}{6 \pi \eta } \left\langle \frac{1}{R_h} \right\rangle $$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 to 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:

$$ D=\frac{D_0}{N} + \frac{k_B T}{6 \pi \eta } \left\langle \frac{1}{R_h} \right\rangle \left( 1- \left\langle\frac{R_h}{L} \right\rangle \right) $$It is quite difficult to fit this analytical expression to simulation data with good accuracy.
It will need a LB fluid, long simulation times and a careful analysis. For this tutorial, we
will use an implicit solvent and short polymer lengths to keep the runtime short. If you want
to collect data suitable for the Zimm model, simply set the global variable `POLYMER_MODEL`

to
`'Zimm'`

, and increase the box size and number of beads in the polymer.

We want to determine the long-time self diffusion coefficient from the mean square displacement of the center-of-mass of a single polymer. For large $t$ the mean square displacement is proportional to the time and the diffusion coefficient occurs as a prefactor:

$$ D = \lim_{t\to\infty}\left[ \frac{1}{6t} \left\langle \left(\vec{r}(t) - \vec{r}(0)\right)^2 \right\rangle \right]. $$This equation can be found in virtually any simulation textbook, like [4]. We will set up a polymer in an implicit solvent, 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 will have a couple of steps inbetween and divide the full problem into subproblems that allow to (hopefully) fully understand the process.

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

```
system.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=1.0, sigma=1.0, shift=0.25, cutoff=1.1225)
```

creates a Lennard-Jones interaction with $\varepsilon=1.$, $\sigma=1.$, $r_{\text{cut}} = 1.1225$ and $\varepsilon_{\text{shift}}=0.25$ between particles of type 0, which is the desired repulsive interaction. The command

```
fene = espressomd.interactions.FeneBond(k=7, r_0=1, d_r_max=2)
```

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

```
system.bonded_inter.add(fene)
```

and to apply the bonded interaction to all monomer pairs of the polymer as shown in the script below.

**ESPResSo** provides a function that tries to find monomer positions that minimize the overlap between
monomers of a chain, *e.g.*:

```
positions = espressomd.polymer.linear_polymer_positions(n_polymers=1,
beads_per_chain=10,
bond_length=1, seed=42,
min_distance=0.9)
```

which would create positions for a single polymer with 10 monomers. Please check the documentation for a more detailed description.

The first task is to compute the average hydrodynamic radius $R_h$, end-to-end distance $R_F$ and radius of gyration $R_g$ for different polymer lengths. This will be achieved with the corresponding observables described in the user guide under Analysis / Direct analysis routines / Chains.

The second task is to estimate the polymer diffusion coefficient for different polymer lengths using two methods:

- the center of mass mean squared displacement method (introduced in a previous part of this tutorial)
- the center of mass velocity autocorrelation method (also known as Green–Kubo method)

For this purpose we can again use the multiple tau correlator.

`build_polymer(system, n_monomers, polymer_params, fene)`

that creates
a linear polymer made of `n_monomers`

particles, with parameters `polymer_params`

. The particles need
to be created and linked together with the `fene`

bond.

In [1]:

```
def build_polymer(system, n_monomers, polymer_params, fene):
positions = espressomd.polymer.linear_polymer_positions(
beads_per_chain=n_monomers, **polymer_params)
p_previous = None
for i, pos in enumerate(positions[0]):
p = system.part.add(pos=pos)
if p_previous is not None:
p.add_bond((fene, p_previous))
p_previous = p
```

`correlator_msd(pids_monomers, tau_max)`

that returns a center-of-mass
mean-squared displacement correlator that is updated every time step, and a function with signature
`correlator_gk(pids_monomers, tau_max)`

that returns a center-of-mass velocity correlator that is updated
every 10 time steps. You can find examples in the user guide section
calculating a particle's diffusion coefficient.

In [2]:

```
def correlator_msd(pids_monomers, tau_max):
com_pos = espressomd.observables.ComPosition(ids=pids_monomers)
com_pos_cor = espressomd.accumulators.Correlator(
obs1=com_pos, tau_lin=16, tau_max=tau_max, delta_N=5,
corr_operation="square_distance_componentwise", compress1="discard1")
return com_pos_cor
def correlator_gk(pids_monomers, tau_max):
com_vel = espressomd.observables.ComVelocity(ids=pids_monomers)
com_vel_cor = espressomd.accumulators.Correlator(
obs1=com_vel, tau_lin=16, tau_max=tau_max, delta_N=10,
corr_operation="scalar_product", compress1="discard1")
return com_vel_cor
```

In [3]:

```
def solvent_langevin(system, kT, gamma):
'''
Implicit solvation model based on Langevin dynamics (Rouse model).
'''
system.thermostat.set_langevin(kT=kT, gamma=gamma, seed=42)
def solvent_lbm(system, kT, gamma):
'''
Lattice-based solvation model based on the LBM (Zimm model).
'''
lbf = espressomd.lb.LBFluidWalberla(kT=kT, seed=42, agrid=1, density=1,
kinematic_viscosity=5, tau=system.time_step,
single_precision=True)
system.lb = lbf
system.thermostat.set_lb(LB_fluid=lbf, gamma=gamma, seed=0)
```

In [4]:

```
import logging
import tqdm
import sys
import numpy as np
import scipy.optimize
import espressomd
import espressomd.analyze
import espressomd.accumulators
import espressomd.observables
import espressomd.polymer
logging.basicConfig(level=logging.INFO, stream=sys.stdout)
espressomd.assert_features(['LENNARD_JONES'])
# Setup constants
BOX_L = 16
TIME_STEP = 0.01
LOOPS = 4000
STEPS = 200
KT = 1.0
GAMMA = 5.0
POLYMER_PARAMS = {'n_polymers': 1, 'bond_length': 1, 'seed': 42, 'min_distance': 0.9}
POLYMER_MODEL = 'Rouse'
assert POLYMER_MODEL in ('Rouse', 'Zimm')
if POLYMER_MODEL == 'Zimm':
espressomd.assert_features(['WALBERLA'])
import espressomd.lb
# System setup
system = espressomd.System(box_l=3 * [BOX_L])
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, r_0=1, d_r_max=2)
system.bonded_inter.add(fene)
N_MONOMERS = [6,8,10,12,14]
com_pos_tau_results = []
com_pos_msd_results = []
com_vel_tau_results = []
com_vel_acf_results = []
rh_results = []
rf_results = []
rg_results = []
for index, N in enumerate(N_MONOMERS):
logging.info(f"Polymer size: {N}")
build_polymer(system, N, POLYMER_PARAMS, fene)
logging.info("Warming up the polymer chain.")
system.time_step = 0.002
system.integrator.set_steepest_descent(
f_max=1.0,
gamma=10,
max_displacement=0.01)
system.integrator.run(2000)
system.integrator.set_vv()
logging.info("Warmup finished.")
logging.info("Equilibration.")
system.time_step = TIME_STEP
system.thermostat.set_langevin(kT=1.0, gamma=50, seed=42)
system.integrator.run(2000)
logging.info("Equilibration finished.")
system.thermostat.turn_off()
if POLYMER_MODEL == 'Rouse':
solvent_langevin(system, KT, GAMMA)
elif POLYMER_MODEL == 'Zimm':
solvent_lbm(system, KT, GAMMA)
logging.info("Warming up the system with the fluid.")
system.integrator.run(1000)
logging.info("Warming up the system with the fluid finished.")
# configure MSD correlator
com_pos_cor = correlator_msd(np.arange(N), LOOPS * STEPS)
system.auto_update_accumulators.add(com_pos_cor)
# configure Green-Kubo correlator
com_vel_cor = correlator_gk(np.arange(N), LOOPS * STEPS)
system.auto_update_accumulators.add(com_vel_cor)
logging.info("Sampling started.")
rhs = np.zeros(LOOPS)
rfs = np.zeros(LOOPS)
rgs = np.zeros(LOOPS)
for i in tqdm.trange(LOOPS):
system.integrator.run(STEPS)
rhs[i] = system.analysis.calc_rh(
chain_start=0,
number_of_chains=1,
chain_length=N)[0]
rfs[i] = system.analysis.calc_re(
chain_start=0,
number_of_chains=1,
chain_length=N)[0]
rgs[i] = system.analysis.calc_rg(
chain_start=0,
number_of_chains=1,
chain_length=N)[0]
logging.info("Sampling finished.")
# store results
com_pos_cor.finalize()
com_pos_tau_results.append(com_pos_cor.lag_times())
com_pos_msd_results.append(np.sum(com_pos_cor.result(), axis=1))
com_vel_cor.finalize()
com_vel_tau_results.append(com_vel_cor.lag_times())
com_vel_acf_results.append(com_vel_cor.result())
rh_results.append(rhs)
rf_results.append(rfs)
rg_results.append(rgs)
# reset system
system.part.clear()
system.lb = None
system.thermostat.turn_off()
system.auto_update_accumulators.clear()
rh_results = np.array(rh_results)
rf_results = np.array(rf_results)
rg_results = np.array(rg_results)
com_pos_tau_results = np.array(com_pos_tau_results)
com_pos_msd_results = np.reshape(com_pos_msd_results, [len(N_MONOMERS), -1])
com_vel_tau_results = np.array(com_vel_tau_results)
com_vel_acf_results = np.reshape(com_vel_acf_results, [len(N_MONOMERS), -1])
```

INFO:root:Polymer size: 6 INFO:root:Warming up the polymer chain. INFO:root:Warmup finished. INFO:root:Equilibration. INFO:root:Equilibration finished. INFO:root:Warming up the system with the fluid. INFO:root:Warming up the system with the fluid finished. INFO:root:Sampling started.

100%|██████████| 4000/4000 [00:20<00:00, 191.83it/s]

INFO:root:Sampling finished. INFO:root:Polymer size: 8 INFO:root:Warming up the polymer chain. INFO:root:Warmup finished. INFO:root:Equilibration.

INFO:root:Equilibration finished. INFO:root:Warming up the system with the fluid. INFO:root:Warming up the system with the fluid finished. INFO:root:Sampling started.

100%|██████████| 4000/4000 [00:21<00:00, 182.20it/s]

INFO:root:Sampling finished. INFO:root:Polymer size: 10 INFO:root:Warming up the polymer chain. INFO:root:Warmup finished. INFO:root:Equilibration.

100%|██████████| 4000/4000 [00:22<00:00, 174.98it/s]

INFO:root:Sampling finished. INFO:root:Polymer size: 12 INFO:root:Warming up the polymer chain. INFO:root:Warmup finished. INFO:root:Equilibration.

100%|██████████| 4000/4000 [00:23<00:00, 169.58it/s]

INFO:root:Sampling finished. INFO:root:Polymer size: 14 INFO:root:Warming up the polymer chain.

INFO:root:Warmup finished. INFO:root:Equilibration. INFO:root:Equilibration finished. INFO:root:Warming up the system with the fluid. INFO:root:Warming up the system with the fluid finished. INFO:root:Sampling started.

100%|██████████| 4000/4000 [00:24<00:00, 160.80it/s]

INFO:root:Sampling finished.

In [5]:

```
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
plt.rcParams.update({'font.size': 18})
```

In [6]:

```
def standard_error_mean_autocorrelation(time_series, variable_label):
'''
Calculate the mean and the correlation-corrected standard error
of the mean of time series by integrating the autocorrelation
function. See Janke 2002 [5] and Weigel, Janke 2010 [6].
Due to the short simulation length, it is not possible to fit an
exponential to the long-time tail. Instead, return a percentile.
'''
summary = []
fig = plt.figure(figsize=(10, 6))
for signal, N in zip(time_series, N_MONOMERS):
acf = espressomd.analyze.autocorrelation(signal - np.mean(signal))
# the acf cannot be integrated beyond tau=N/2
integral = np.array([acf[0] + 2 * np.sum(acf[1:j]) for j in np.arange(1, len(acf) // 2)])
# remove the noisy part of the integral
negative_number_list = np.nonzero(integral < 0)
if negative_number_list[0].size:
integral = integral[:int(0.95 * negative_number_list[0][0])]
# compute the standard error of the mean
std_err = np.sqrt(integral / acf.size)
# due to the small sample size, the long-time tail is not
# well resolved and cannot be fitted, so we use a percentile
asymptote = np.percentile(std_err, 75)
# plot the integral and asymptote
p = plt.plot([0, len(std_err)], 2 * [asymptote], '--')
plt.plot(np.arange(len(std_err)) + 1, std_err,
'-', color=p[0].get_color(),
label=rf'$\int {variable_label}$ for N={N}')
summary.append((np.mean(signal), asymptote))
plt.xlabel(r'Lag time $\tau / \Delta t$')
plt.ylabel(rf'$\int_{{-\tau}}^{{+\tau}} {variable_label}$')
plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 0))
plt.legend()
plt.show()
return np.array(summary)
def fitting_polymer_theory(polymer_model, n_monomers, diffusion, rh_exponent):
'''
Fit the appropriate polymer diffusion coefficient equation (Rouse or
Kirkwood-Zimm).
'''
def rouse(x, a):
return a / x
def kirkwood_zimm(x, a, b, exponent):
return a / x + b / x**exponent
x = np.linspace(min(n_monomers) - 0.5, max(n_monomers) + 0.5, 20)
if polymer_model == 'Rouse':
popt, _ = scipy.optimize.curve_fit(rouse, n_monomers, diffusion)
label = rf'$D^{{\mathrm{{fit}}}} = \frac{{{popt[0]:.2f}}}{{N}}$'
y = rouse(x, popt[0])
elif polymer_model == 'Zimm':
fitting_function = kirkwood_zimm
popt, _ = scipy.optimize.curve_fit(
lambda x, a, b: kirkwood_zimm(x, a, b, rh_exponent), n_monomers, diffusion)
y = kirkwood_zimm(x, popt[0], popt[1], rh_exponent)
label = f'''\
$D^{{\\mathrm{{fit}}}} = \
\\frac{{{popt[0]:.2f}}}{{N}} + \
\\frac{{{popt[1] * 6 * np.pi:.3f} }}{{6\\pi}} \\cdot \
\\frac{{{1}}}{{N^{{{rh_exponent:.2f}}}}}$ \
'''
return x, y, label, popt
```

How do $R_h$, $R_g$, $R_F$ and the diffusion coefficient $D$ depend on the number of monomers? You can refer to the Flory theory of polymers, and assume you are simulating a real polymer in a good solvent, with Flory exponent $\nu \approx 0.588$.

In [7]:

```
rf_summary = standard_error_mean_autocorrelation(rf_results, r'\operatorname{acf}(R_F)')
rf_exponent, rf_prefactor = np.polyfit(np.log(N_MONOMERS), np.log(rf_summary[:, 0]), 1)
rf_prefactor = np.exp(rf_prefactor)
fig = plt.figure(figsize=(10, 6))
x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)
plt.plot(x, rf_prefactor * x**rf_exponent, '-',
label=rf'$R_F^{{\mathrm{{fit}}}} = {rf_prefactor:.2f} N^{{{rf_exponent:.2f}}}$')
plt.errorbar(N_MONOMERS, rf_summary[:, 0],
yerr=rf_summary[:, 1],
ls='', marker='o', capsize=5, capthick=1,
label=r'$R_F^{\mathrm{simulation}}$')
plt.xlabel('Number of monomers $N$')
plt.ylabel(r'End-to-end distance [$\sigma$]')
plt.legend()
plt.show()
```

In [8]:

```
rg_summary = standard_error_mean_autocorrelation(rg_results, r'\operatorname{acf}(R_g)')
rg_exponent, rg_prefactor = np.polyfit(np.log(N_MONOMERS), np.log(rg_summary[:, 0]), 1)
rg_prefactor = np.exp(rg_prefactor)
fig = plt.figure(figsize=(10, 6))
x = np.linspace(min(N_MONOMERS) - 0.5, max(N_MONOMERS) + 0.5, 20)
plt.plot(x, rg_prefactor * x**rg_exponent, '-',
label=rf'$R_g^{{\mathrm{{fit}}}} = {rg_prefactor:.2f} N^{{{rg_exponent:.2f}}}$')
plt.errorbar(N_MONOMERS, rg_summary[:, 0],
yerr=rg_summary[:, 1],
ls='', marker='o', capsize=5, capthick=1,
label=r'$R_g^{\mathrm{simulation}}$')
plt.xlabel('Number of monomers $N$')
plt.ylabel('Radius of gyration [$\sigma$]')
plt.legend()
plt.show()
```