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, cutof=1.226)
creates a Lennard-Jones interaction with $\varepsilon=1.$, $\sigma=1.$, $r_\text{cut} = 1.125$ 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, 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=5642,
min_distance=0.9)
which would create positions for a single polymer with 10 monomers. Please check the documentation for a more detailed description.
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.
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.
In [ ]:
import logging
import sys
import numpy as np
import espressomd
import espressomd.accumulators
import espressomd.observables
import espressomd.polymer
import espressomd.minimize_energy
logging.basicConfig(level=logging.INFO, stream=sys.stdout)
espressomd.assert_features(['LENNARD_JONES'])
# Setup constant
TIME_STEP = 0.01
LOOPS = 100
STEPS = 100
# System setup
system = espressomd.System(box_l=[32.0, 32.0, 32.0])
system.cell_system.skin = 0.4
N_MONOMERS = 20 # The number of monomers has been set to be 20 as default
# Change this value for further simulations
# 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)
# Setup polymer of part_id 0 with fene bond
positions = espressomd.polymer.linear_polymer_positions(n_polymers=1,
beads_per_chain=N_MONOMERS,
bond_length=1, seed=5642,
min_distance=0.9)
for i, pos in enumerate(positions[0]):
id = len(system.part)
system.part.add(id=id, pos=pos)
if i > 0:
system.part[id].add_bond((fene, id - 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("LB fluid warming finished.")
# configure correlator
com_pos = espressomd.observables.ComPosition(ids=range(N_MONOMERS))
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.")
for i in range(LOOPS):
system.integrator.run(STEPS)
logging.info("Sampling finished.")
correlator.finalize()
corrdata = correlator.result()
tau = corrdata[:, 0]
msd = corrdata[:, 2] + corrdata[:, 3] + corrdata[:, 4]
rh = system.analysis.calc_rh(chain_start=0, number_of_chains=1, chain_length=N_MONOMERS)[0]
[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 to 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.