In this tutorial we are going to investigate diffusion of a dissolved polymer using 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}} = 2^{1/6} \simeq 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:
For this purpose we can again use the multiple tau correlator.
Write a function with signature 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.
# SOLUTION CELL
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
Write a function with signature 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.
# SOLUTION CELL
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
You can simulate a polymer in the Rouse regime using an implicit solvent model, e.g. Langevin dynamics, or in the Zimm regime using a lattice-Boltzmann fluid.
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)
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
import espressomd.zn
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
MASS = 1.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.time_step = TIME_STEP
system.cell_system.skin = 0.4
# Lennard-Jones interaction
LJ_SIGMA=1.0
LJ_EPSILON=1.0
LJ_CUTOFF=2.0**(1.0 / 6.0)
system.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=LJ_EPSILON, sigma=LJ_SIGMA, shift="auto", cutoff=LJ_CUTOFF)
# Fene interaction
fene = espressomd.interactions.FeneBond(k=7, r_0=1, d_r_max=2)
system.bonded_inter.add(fene)
0
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 = []
### System visualization
print("Simulation video for polymer of size :", N_MONOMERS[2])
# Visualizer Parameters
color = {0: "blue"} # Particle color by type
radii = {0: LJ_SIGMA/2.0} # Particle size by type
# Initializing Visualizer
vis = espressomd.zn.Visualizer(system, colors=color, radii=radii, bonds=True)
Simulation video for polymer of size : 10
for N in N_MONOMERS:
logging.info(f"Polymer size: {N}")
build_polymer(system, N, POLYMER_PARAMS, fene)
logging.info("Removing overlaps ...")
FMAX = 0.001 * LJ_SIGMA * MASS / system.time_step**2
system.integrator.set_steepest_descent(
f_max=FMAX,
gamma=10,
max_displacement=0.01)
system.integrator.run(100)
assert np.all(np.abs(system.part.all().f) < FMAX), "Overlap removal did not converge!"
system.integrator.set_vv()
logging.info("Remove overlap 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)
if N == N_MONOMERS[2]:
vis.update()
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:Removing overlaps ... INFO:root:Remove overlap 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:10<00:00, 389.83it/s]
INFO:root:Sampling finished. INFO:root:Polymer size: 8 INFO:root:Removing overlaps ... INFO:root:Remove overlap 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:11<00:00, 353.75it/s]
INFO:root:Sampling finished. INFO:root:Polymer size: 10 INFO:root:Removing overlaps ... INFO:root:Remove overlap 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:11<00:00, 346.44it/s]
INFO:root:Sampling finished. INFO:root:Polymer size: 12 INFO:root:Removing overlaps ... INFO:root:Remove overlap 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:11<00:00, 342.82it/s]
INFO:root:Sampling finished. INFO:root:Polymer size: 14 INFO:root:Removing overlaps ... INFO:root:Remove overlap 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:12<00:00, 317.59it/s]
INFO:root:Sampling finished.
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
plt.rcParams.update({'font.size': 18})
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$.
Plot the end-to-end distance $R_F$ of the polymer as a function of the number of monomers. What relation do you observe?
The end-to-end distance follows the law $R_F = c_F N^\nu$ with $c_F$ a constant and $\nu$ the Flory exponent.
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()
Plot the radius of gyration $R_g$ of the polymer as a function of the number of monomers. What relation do you observe?
The radius of gyration follows the law $R_g = c_g N^\nu$ with $c_g$ a constant and $\nu$ the Flory exponent.
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(r'Radius of gyration [$\sigma$]')
plt.legend()
plt.show()